Note
Go to the end to download the full example code.
Applying a Filter to a Time-Series#
This example demonstrates low pass filtering a time-series by applying a weighted running mean over the time dimension.
The time-series used is the Darwin-only Southern Oscillation index (SOI), which is filtered using two different Lanczos filters, one to filter out time-scales of less than two years and one to filter out time-scales of less than 7 years.
References#
Duchon C. E. (1979) Lanczos Filtering in One and Two Dimensions. Journal of Applied Meteorology, Vol 18, pp 1016-1022.
Trenberth K. E. (1984) Signal Versus Noise in the Southern Oscillation. Monthly Weather Review, Vol 112, pp 326-332
import matplotlib.pyplot as plt
import numpy as np
import iris
import iris.plot as iplt
def low_pass_weights(window, cutoff):
"""Calculate weights for a low pass Lanczos filter.
Parameters
----------
window : int
The length of the filter window.
cutoff : float
The cutoff frequency in inverse time steps.
"""
order = ((window - 1) // 2) + 1
nwts = 2 * order + 1
w = np.zeros([nwts])
n = nwts // 2
w[n] = 2 * cutoff
k = np.arange(1.0, n)
sigma = np.sin(np.pi * k / n) * n / (np.pi * k)
firstfactor = np.sin(2.0 * np.pi * cutoff * k) / (np.pi * k)
w[n - 1 : 0 : -1] = firstfactor * sigma
w[n + 1 : -1] = firstfactor * sigma
return w[1:-1]
def main():
# Load the monthly-valued Southern Oscillation Index (SOI) time-series.
fname = iris.sample_data_path("SOI_Darwin.nc")
soi = iris.load_cube(fname)
# Window length for filters.
window = 121
# Construct 2-year (24-month) and 7-year (84-month) low pass filters
# for the SOI data which is monthly.
wgts24 = low_pass_weights(window, 1.0 / 24.0)
wgts84 = low_pass_weights(window, 1.0 / 84.0)
# Apply each filter using the rolling_window method used with the weights
# keyword argument. A weighted sum is required because the magnitude of
# the weights are just as important as their relative sizes.
soi24 = soi.rolling_window("time", iris.analysis.SUM, len(wgts24), weights=wgts24)
soi84 = soi.rolling_window("time", iris.analysis.SUM, len(wgts84), weights=wgts84)
# Plot the SOI time series and both filtered versions.
plt.figure(figsize=(9, 4))
iplt.plot(
soi,
color="0.7",
linewidth=1.0,
linestyle="-",
alpha=1.0,
label="no filter",
)
iplt.plot(
soi24,
color="b",
linewidth=2.0,
linestyle="-",
alpha=0.7,
label="2-year filter",
)
iplt.plot(
soi84,
color="r",
linewidth=2.0,
linestyle="-",
alpha=0.7,
label="7-year filter",
)
plt.ylim([-4, 4])
plt.title("Southern Oscillation Index (Darwin Only)")
plt.xlabel("Time")
plt.ylabel("SOI")
plt.legend(fontsize=10)
iplt.show()
if __name__ == "__main__":
main()
Total running time of the script: (0 minutes 0.704 seconds)