Filtering seasonal time series outliers

Question

I have some weather-related data that is highly seasonal, what i want to do is to identify some 'outliers' and change these outliers to reasonable values (i dont wan't to remove them).

I'm using the hampel filter from the pracma package to identify and correct these outliers, but the problem is that some seasonal-peaks are being lowered too much. I guess this function doesn't take into account seasonality. The plot below shows this (window of 24 months), the red line is the filtered data.

Is there a way to correct this?In other words, take the 'peaks' and average them between seasons? Thanks!

Answer 1

As suggested by @Michael, here's some code that remove the seasonal component, apply a outlier filter on the remaining serie, and adds the seasonal component again (use by your own risk):

filter.outliers = function(serie.orig)
{
require(forecast)
require(pracma)

## get a stl decomposition object
stl.decomp = stl(serie.orig, s.window = 'periodic', t.window = 13, robust = T)

## remove the seasonal component
serie.desaz = seasadj(stl.decomp)

## apply a hampel filter
serie.desaz2 = hampel(serie.desaz, k = 12)$y

## put back the seasonal component
return(serie.desaz2 + stl.decomp$time.series[, 'seasonal'])
}