I´d like to extract years from a time series index (the underlying time series is of monthly frequency). The reason I want to do it is creating a yearly axis.

I am using a package Rbeast.

Here is my result for tsp(NDVI_Chobe001.ts) [1] 2002.500 2020.417 12.000

##Here is a more detailed time series

 time(NDVI_Chobe001.ts)
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep
2002                                                       2002.500 2002.583 2002.667
2003 2003.000 2003.083 2003.167 2003.250 2003.333 2003.417 2003.500 2003.583 2003.667
2004 2004.000 2004.083 2004.167 2004.250 2004.333 2004.417 2004.500 2004.583 2004.667
2005 2005.000 2005.083 2005.167 2005.250 2005.333 2005.417 2005.500 2005.583 2005.667
2006 2006.000 2006.083 2006.167 2006.250 2006.333 2006.417 2006.500 2006.583 2006.667
2007 2007.000 2007.083 2007.167 2007.250 2007.333 2007.417 2007.500 2007.583 2007.667
2008 2008.000 2008.083 2008.167 2008.250 2008.333 2008.417 2008.500 2008.583 2008.667
2009 2009.000 2009.083 2009.167 2009.250 2009.333 2009.417 2009.500 2009.583 2009.667
2010 2010.000 2010.083 2010.167 2010.250 2010.333 2010.417 2010.500 2010.583 2010.667
2011 2011.000 2011.083 2011.167 2011.250 2011.333 2011.417 2011.500 2011.583 2011.667
2012 2012.000 2012.083 2012.167 2012.250 2012.333 2012.417 2012.500 2012.583 2012.667
2013 2013.000 2013.083 2013.167 2013.250 2013.333 2013.417 2013.500 2013.583 2013.667
2014 2014.000 2014.083 2014.167 2014.250 2014.333 2014.417 2014.500 2014.583 2014.667
2015 2015.000 2015.083 2015.167 2015.250 2015.333 2015.417 2015.500 2015.583 2015.667
2016 2016.000 2016.083 2016.167 2016.250 2016.333 2016.417 2016.500 2016.583 2016.667
2017 2017.000 2017.083 2017.167 2017.250 2017.333 2017.417 2017.500 2017.583 2017.667
2018 2018.000 2018.083 2018.167 2018.250 2018.333 2018.417 2018.500 2018.583 2018.667
2019 2019.000 2019.083 2019.167 2019.250 2019.333 2019.417 2019.500 2019.583 2019.667
2020 2020.000 2020.083 2020.167 2020.250 2020.333 2020.417                           
          Oct      Nov      Dec
2002 2002.750 2002.833 2002.917
2003 2003.750 2003.833 2003.917
2004 2004.750 2004.833 2004.917
2005 2005.750 2005.833 2005.917
2006 2006.750 2006.833 2006.917
2007 2007.750 2007.833 2007.917
2008 2008.750 2008.833 2008.917
2009 2009.750 2009.833 2009.917
2010 2010.750 2010.833 2010.917
2011 2011.750 2011.833 2011.917
2012 2012.750 2012.833 2012.917
2013 2013.750 2013.833 2013.917
2014 2014.750 2014.833 2014.917
2015 2015.750 2015.833 2015.917
2016 2016.750 2016.833 2016.917
2017 2017.750 2017.833 2017.917
2018 2018.750 2018.833 2018.917
2019 2019.750 2019.833 2019.917
2020                     

##Here is the dput results Also this is the result of dput of the time series

dput(NDVI_Chobe001.ts)
structure(c(0.258672185, 0.237639041, 0.223988035, 0.22397988, 
0.28132144, 0.387710909, 0.556186453, 0.719580311, 0.443294248, 
0.314433357, 0.292755672, 0.278023297, 0.246809774, 0.230477039, 
0.228955071, 0.234966762, 0.396659718, 0.491330645, 0.716274496, 
0.73417416, 0.576279481, 0.479403466, 0.423930923, 0.377550602, 
0.313801774, 0.297335261, 0.27886131, 0.285054814, 0.323942137, 
0.482912961, 0.544500134, 0.692623308, 0.512637643, 0.526592884, 
0.432898755, 0.331892323, 0.294543398, 0.274633904, 0.247602217, 
0.24889248, 0.268682448, 0.623132479, 0.732958587, 0.794572999, 
0.697092229, 0.588587711, 0.389213056, 0.348485514, 0.282932917, 
0.264831569, 0.254937301, 0.277647575, 0.342873991, 0.450392312, 
0.631838679, 0.716638566, 0.508704542, 0.423229761, 0.324663042, 
0.321186864, 0.295935102, 0.255618784, 0.257619908, 0.27254904, 
0.270580191, 0.432220414, 0.585539863, 0.716581759, 0.593955107, 
0.559614127, 0.388112159, 0.386833323, 0.323056858, 0.296201373, 
0.299547175, 0.340499135, 0.487932489, 0.531619232, 0.755898976, 
0.630418376, 0.624101053, 0.46624194, 0.407831925, 0.396696504, 
0.348860597, 0.311130303, 0.301932775, 0.359536613, 0.389388896, 
0.639178419, 0.693073002, 0.544238686, 0.608043068, 0.520223438, 
0.489590537, 0.371915235, 0.322345492, 0.285747424, 0.262060636, 
0.290601893, 0.272739968, 0.465184219, 0.597999142, 0.58280379, 
0.498312536, 0.351555151, 0.313456794, 0.30176279, 0.272389062, 
0.256200802, 0.257570355, 0.261360949, 0.280664516, 0.472871998, 
0.635346196, 0.809469884, 0.708410897, 0.454619991, 0.374153554, 
0.31216354, 0.277643235, 0.273912332, 0.279202729, 0.274954368, 
0.3220379, 0.66542086, 0.786922135, 0.749038686, 0.494324261, 
0.380828443, 0.31699487, 0.293321759, 0.275010923, 0.263032267, 
0.254334753, 0.270928799, 0.319622211, 0.699834174, 0.733615599, 
0.743950233, 0.701439492, 0.514223884, 0.414668945, 0.353937354, 
0.304874948, 0.263697644, 0.25397833, 0.270116643, 0.310931558, 
0.529621168, 0.769325391, 0.676328487, 0.561630209, 0.550257129, 
0.408412794, 0.362071348, 0.294350827, 0.271100888, 0.26194204, 
0.265314367, 0.311826075, 0.556399621, 0.586507353, 0.715697274, 
0.71326184, 0.554275685, 0.423167172, 0.332214247, 0.316422341, 
0.264796933, 0.253106015, 0.274547234, 0.274233387, 0.331729551, 
0.568674796, 0.678751599, 0.651238604, 0.550351642, 0.50121297, 
0.39598441, 0.373748125, 0.342521719, 0.308862893, 0.368914514, 
0.346751371, 0.590197072, 0.550842239, 0.641805925, 0.73276961, 
0.616976479, 0.501595155, 0.441702349, 0.420716154, 0.297677153, 
0.302351869, 0.347572841, 0.347128221, 0.488411226, 0.739632011, 
0.773688839, 0.527578039, 0.531276077, 0.481584383, 0.450086834, 
0.331415825, 0.298545112, 0.281891087, 0.301013691, 0.334150401, 
0.537372378, 0.756594273, 0.778707894, 0.728867792, 0.634829201, 
0.415475576, 0.353712963), .Tsp = c(2002.5, 2020.41666666667, 
12), class = "ts")

Here is my code

#Here I am Specifying the option parameters explicitly opt=list() #Create an empty list to append individual model parameters

    opt$period=12           #Period of the cyclic/seasonal component of the modis time series
    opt$startTime=2002.500000   ##specify start time
    
    

Here I am plotting the data

    out=beast(NDVI_Chobe001.ts, opt)
    p<-plot(out)
    

However, the plot time series is in monthly frequency , see example here[enter image description here][1]

##Now I am extracting years from the time series to use when plotting

    time <- as.yearmon( time(NDVI_Chobe001.ts))
    time<-format.Date(x,"%Y")
    
    plot(out,axis(1,time, labels = TRUE))
    

But this is not working, it is still displayed in monthly frequency, How do I lot my time series by year?

      [1]: https://i.stack.imgur.com/VNxJC.jpg

Solution 1:

I know this question is old and my late reply may be irrelevant. Regardless, in the old version of Rbeast, if given a ts input, the beast function uses only its data vector not the time tag information (i.e., .tsp attributes such as start and frequency). In other words, the time info has to be additionally provided to beast; if not, the indices 1:N are used as default times. In the latest version (Rbeast v0.9.2), I changed the API a little bit; now beast can handle the time tag of ts objects.

Below is a test using your sample data for time series decomposition and changepoint detection.

ndvi= c( 0.258,0.237,0.223,0.223,0.281,0.387,0.556,0.719,0.443,0.314,0.292,0.278,0.246,0.230,0.228,0.234,0.396,0.491,0.716, 0.734,0.576,0.479,0.423,0.377,0.313,0.297,0.278,0.285,0.323,0.482,0.544,0.692,0.512,0.526,0.432,0.331,0.294,0.274,
     0.247,0.248,0.268,0.623,0.732,0.794,0.697,0.588,0.389,0.348,0.282,0.264,0.254,0.277,0.342,0.450,0.631,0.716,0.508, 0.423,0.324,0.321,0.295,0.255,0.257,0.272,0.270,0.432,0.585,0.716,0.593,0.559,0.388,0.386,0.323,0.296,0.299,0.340,
     0.487,0.531,0.755,0.630,0.624,0.466,0.407,0.396,0.348,0.311,0.301,0.359,0.389,0.639,0.693,0.544,0.608,0.520,0.489, 0.371,0.322,0.285,0.262,0.290,0.272,0.465,0.597,0.582,0.498,0.351,0.313,0.301,0.272,0.256,0.257,0.261,0.280,0.472,
     0.635,0.809,0.708,0.454,0.374,0.312,0.277,0.273,0.279,0.274,0.322,0.665,0.786,0.749,0.494,0.380,0.316,0.293,0.275, 0.263,0.254,0.270,0.319,0.699,0.733,0.743,0.701,0.514,0.414,0.353,0.304,0.263,0.253,0.270,0.310,0.529,0.769,0.676,
     0.561,0.550,0.408,0.362,0.294,0.271,0.261,0.265,0.311,0.556,0.586,0.715,0.713,0.554,0.423,0.332,0.316,0.264,0.253, 0.274,0.274,0.331,0.568,0.678,0.651,0.550,0.501,0.395,0.373,0.342,0.308,0.368,0.346,0.590,0.550,0.641,0.732,0.616,
     0.501,0.441,0.420,0.297,0.302,0.347,0.347,0.488,0.739,0.773,0.527,0.531,0.481,0.450,0.331,0.298,0.281,0.301,0.334, 0.537,0.756,0.778,0.728,0.634,0.415,0.353)

ndvi is a data vector; create a ts object out of it as the input to beast:

ndvi_ts = ts(ndvi, start=2002.5, frequency = 12)

out = beast(ndvi_ts) 
plot(out)

# By default, seasonality is fitted as a harmonic curve; use SVD-based bases instead
out = beast(ndvi_ts, season='svd')
plot(out)

# If the input is a pure data vector, the time info has to be specified 
# via 'start', 'deltat' (delta T), and 'freq'; otherwise, default times 
# (i.e., 1:216) would be used.
out = beast( ndvi, start=2002.5, deltat=1/12, freq=12, season='svd')
plot(out)

The beast output is a LIST object of class 'beast' and can't be directly plotted with ggplot2. The plot function above is not R's own function but Rbeast's implementation of the S3 method 'plot.beast'. If customary plots are needed, individual variables of the output can be extracted for use in R's base plot or ggplot2. Below is an example:

out=beast(ndvi_ts, season='svd')

# out is a LIST variable; extract some elements
t    = out$time            # time
y    = out$trend$Y         # fitted trend compnt
yci  = out$trend$CI[,,1]   # estimated 95% CI interval for trend
cp   = out$trend$cp        # most probable changepoints in trend 
ncp  = sum(!is.nan(cp))    # number of changepoints given in cp

ymax=max(yci)
ymin=min(yci)    
plot( x=c(min(t),max(t)), y=c(ymin,ymax), type='n', xlab='time', ylab='trend'  )
lines(t,y)
polygon( x=c(t,rev(t)), y=c( yci[,1], rev(yci[,2]) ) , col=rgb(0,1,0,.2),border=NA )

yrange=par('usr')[3:4]
for (i in 1:ncp) {
  lines( c(cp[i], cp[i] ), y=yrange, type='l', lty='dashed' )  
}

If multiple subplots are needed, one way is to use R's par(new=TRUE) option, as shown below.

t    = out$time            # time
y    = out$trend$Y         # fitted trend compnt
yci  = out$trend$CI[,,1]   # estimated 95% CI interval for trend
cp   = out$trend$cp        # most probable changepoints in trend 
ncp  = sum(!is.nan(cp))    # number of changepoints given in cp
prob = out$trend$cpOccPr   # probability of changepoints occurrence over time

ymax=max(yci)
ymin=min(yci)

plot.new()

# FIRST PLOT: Fine-tune the four margin numbers in plt to re-position
par(plt = c(0.15, 1-0.01, 1-0.01-.5, 1-0.01-.5+0.4),new=TRUE) 

plot( x=c(min(t),max(t)), y=c(ymin,ymax), type='n', xaxt='n', xlab=NA,  ylab='trend')
lines( t, y )
polygon( x=c( t, rev(t) ), y=c( yci[,1], rev(yci[,2]) ) , col=rgb(0,1,0,.2),border=NA )

yrange=par('usr')[3:4]
for (i in 1:ncp) {
  lines( c(cp[i], cp[i] ), y=yrange, type='l', lty='dashed' )  
}

# SECOND PLOT: Fine-tune the four margin numbers in plt to re-position
par(plt = c(0.15, 1-0.01, 1-0.01-.5-0.3, 1-0.01-.5-0.3+0.3),new=TRUE)
plot(  x=t, y=prob, type='l', xlab='time', ylab='Prob'  )
     

Not sure how useful this is. In any case, thanks for testing out Rbeast. If new features are desired, pls let me know.