plot.bfast: Better plotting

[GreatEmerald]

The current plot method in plot.bfast shows the components one by one, but there is no way to plot everything in one plot. It would be very convenient to have that.

Here's a small snippet of my custom code that produces pretty nice results:

MyTS = bfastts(ChangedEVITS[9,], dates, "10-day") # Some time series
MyTS = na.approx(MyTS) # This step is optional, just for break-free visualisation
MyBF = bfast(MyTS, GetBreakNumber(dates), season="harmonic", max.iter=3)
MyIter = length(MyBF$output)

MyFit = MyBF$output[[MyIter]]$St + MyBF$output[[MyIter]]$Tt
MyTrend = MyBF$output[[MyIter]]$Tt
MyTimes = as.numeric(time(MyBF$output[[MyIter]]$Tt))
MyTimes = MyTimes[!is.na(MyBF$output[[MyIter]]$Tt)]
MyBreaks = MyTimes[MyBF$output[[MyIter]]$bp.Vt$breakpoints]
plot(MyTS); lines(MyFit, col="blue"); lines(MyTrend, col="green"); abline(v=MyBreaks, col="red")

Sample output:

What's missing is a legend. Suggestions are also welcome!