season estimates anomalies associated with the annual cycle or of a time series.
[As,ts] = season(A,t) [As,ts] = season(...,'daily') [As,ts] = season(...,'monthly') [As,ts] = season(...,'detrend',DetrendOption) [As,ts] = season(...,'dim',dimension) As = season(...,'full')
[As,ts] = season(A,t) gives the typical seasonal (aka annual) cycle of the time series A corresponding to times t in datenum or datetime format. If t is daily, outputs ts is 1 to 366 and As will contain average values for each of the 366 days of the year. If inputs are monthly, ts is 1:12 and As will contain average values for each of the 12 months of the year.
[As,ts] = season(...,'daily') specifies directly that inputs are daily resolution. The season function will typically figure this out automatically, but if you have large missing gaps in your data you may wish to ensure correct results by specifying daily.
[As,ts] = season(...,'monthly') as above, but forces monthly solution.
[As,ts] = season(...,'detrend',DetrendOption) specifies a baseline relative to which seasonal anomalies are determined. Options are 'linear', 'quadratic', or 'none'. By default, anomalies are calculated after removing the linear least squares trend, but if, for example, warming is strongly nonlinear, you may prefer the 'quadratic' option. Default is 'linear'.
[As,ts] = season(...,'dim',dimension) specifies a dimension along which to assess seasons. By default, if A is 1D, seasonal cycle is returned along the nonsingleton dimension; if A is 2D, season is performed along dimension 1 (time marches down the rows); if A is 3D, season is performed along dimension 3.
As = season(...,'full') returns As for the entire time series A. This is a convenient option for looking at the components of a long time series separately.
Example 1: A single time series
Consider this time series. Let's call it temperature:
t = datenum('jan 1, 1979'):datenum('dec 31, 2016'); % Sinusoidal seasonal signal + warming trend + noise + some mean value: T = 4*sin(doy(t,'decimalyear')*2*pi) + (t-min(t))/5e3 + randn(size(t)) + 15; plot(t,T) axis tight box off datetick('x','keeplimits')
You can see the temperature has a seasonal cycle, but there's a gradual warming over time and there's also some random noise in there. To see the structure of the seasonal cycle a little better, we can look at the same time series as a function of day of year using the doy function. We'll also use the cbdate function to format the colorbar dates.
figure scatter(doy(t),T,10,t,'filled') axis tight colorbar cbdate('yyyy') xlabel 'day of year' ylabel 'temperature' cmocean amp
And again you can see the temperature rising over time, but by plotting as a function of doy we get to see the structure of the seasonal cycle a little better. Now what if we want to separate the effects of long-term warming from the seasonal cycle?
Ts = season(T,t); figure plot(1:366,Ts) box off axis tight xlabel 'day of year' ylabel 'seasonal temperature cycle'
Above we see the seasonal cycle that has an amplitude of 4 degrees, just as we defined it. Some noise remains because each day of the year represents the average of 38 noisy measurements. Nonetheless, it's our best approximation of the seasonal cycle, and it's pretty close.
Example 2: Multiple time series at once
Let's say you have several rows of data, each containing a different time series. The data might look like this (we're building on the time series T from above)
T2 = [T + (t-min(t))/2e3; T - 15*cos(doy(t,'decimalyear')*2*pi) + (t-min(t))/1e2; T + ((t-min(t))/2.5e3).^3.1]; figure plot(t,T2) axis tight box off datetick('x','keeplimits') legend('first row of T2','second row of T2','third row of T2','location','northwest') legend boxoff
What's the annual cycle like for each of those time series? Use season, but specify dimension 2 because each row of T2 contains a time series. If each column of T2 contained a time series, the default dimension 1 would be what we want.
Let's also get the whole time series with the 'full' option rather than the default days 1 to 366.
T2_full = season(T2,t,'dimension',2,'full'); figure plot(t,T2_full) axis([datenum('jan 1, 1995') datenum('jan 1, 2010') -16 16]) box off datetick('x','keeplimits') legend('first row of T2','second row of T2','third row of T2','location','northwest')
Above we see that season has extracted the seasonal cycles of each time series in T2. Even row 3 of T2, which has a long term trend that scales to the power of 3.1 is sufficiently analyzed using the default linear detrending. However, if linear detrending is insufficient for your application, use the 'detrend','quadratic' option to step it up a notch.
Removing seasonal cycles
For convenience, the Climate Data Toolbox has a deseason function, but it's worth noting here, that the deseason function removes seasonal cycles just by subtracting the seasonal component from the original signals, like this:
T_deseasoned = T2 - T2_full; figure plot(t,T2) hold on plot(t,T_deseasoned) axis tight box off datetick('x','keeplimits') legend('first row of T2','second row of T2','third row of T2','location','northwest') legend boxoff
Example 3: Real Sea Ice Data
The examples above used some artifical data we created. But how does this function apply to real data? Let's compare the northern and southern hemisphere seasonal cycles of sea ice extent. Here's what the full time series look like:
load seaice_extent figure plot(t,extent_N,'b') hold on plot(t,extent_S,'r') axis tight box off legend('north','south','location','northwest')
So what's a typical year of sea ice look like?
figure plot(season(extent_N,t),'b') hold on plot(season(extent_S,t),'r') axis tight box off xlabel 'day of year' ylabel('sea ice extent anomaly (10^6 km^2)') ntitle(' typical seasonal cycle ') legend('north','south','location','northwest')
Example 4: Gridded reanalysis data
Let's explore the typical seasonal cycles of sea surface temperatures. Here's a monthly dataset of dimensions 60x55x802:
load pacific_sst sst_season = season(sst,t);
That gives us sst_season, which is 60x55x12, representing one 60x55 grid for each month of the year. Take the difference between the typical maximum and minimum temperatures to get a sense of how strong the seasons are in different parts of the Pacific:
T_range = max(sst_season,,3)-min(sst_season,,3); figure imagescn(lon,lat,T_range) cb = colorbar; ylabel(cb,'magnitude of seasonal cycle (\circ C)') cmocean amp
That's not surprising--Seasonality is strongest close to the poles and in shallow waters.
How this function works
By default, the season function removes the linear trend from a time series, then determines the climatology of the time series by averaging all detrended data for each day (or month) of the year.
For daily data, Leap Years present a challenge. This function calculates the average temperature for February 28th by averaging all data corresponding to the 59th day of the year. But what to do about March 1st? Most years March 1st is the 60th day of the year, but every fourth year March 1st is the 61st day of the year.
This function converts all dates to the day of the year, and calculates seasonal cycles as a function of day of the year. Since there are far fewer years with 366 days, applying the same practice of averaging all data from day 366 would make day 366 sensitive to outliers. Thus, for day 366, the season function averages all data from day 365, day 366, and day 1 of all years.
Seasons vs Climatology
CDT has a function called season and another function called climatology. The only difference is that the output of climatology includes the mean of the variable, whereas the output of season will always have a 0 mean value. Accordingly, the deseason function removes the seasonal component of variability while preserving the overall mean and trends.
In general, CDT assumes that a multi-year record of a variable sampled at subannual resolution can described by
y = y_0 + y_tr + y_season + y_var + y_noise
- y_0 is the long-term mean,
- y_trend is the long-term trend,
- y_season is the typical seasonal anomaly, which the season function obtains after detrending and removing the mean,
- y_var represents interannual variability, and
- y_noise is everything else
In this model,
y_climatology = y_0 + y_season
Other ways to define seasonality
This function was written by Chad A. Greene of the University of Texas Institute for Geophysics (UTIG), July 2017.