R/plot.SPTMres.R

In ick003/SpTMixture: Mixture models with spatio-temporal dependencies

#' Update the seasons for LU
#'
#' Update the seasons for LU
#' @param season.sites the season calculated for each site.
#' @param p.sites the composition of each site.
#' @keywords updating season
#' @export
#' @examples
#' data("wq_analysis_week2")
#' SPTMData(wq.raw.obs, frequency = "quarter")
plot.sptmRes = function(SPTMresobj){

  theta = SPTMresobj$Result$theta
  p.sites = SPTMresobj$Result$latent
  seasonSites = SPTMresobj$seasonSites
  trendSites = SPTMresobj$trendSites

  df = SPTMresobj$data

  y = df$obs.data

  y.p = meanFunc(theta, n.mixt=5, p.sites=p.sites, season.sites = seasonSites, trend.sites = trendSites)

  SiteID = unique(df$obs.data$ID)
  par(mfrow=c(4,4), mar=c(2,2,2,1))
  for(i in 1:16){
    obs = y$obs[y$ID %in% SiteID[i]]
    idx.k = which.max(p.sites[i,])
    pred= y.p[,idx.k]
    dates = y$date[y$ID %in% SiteID[i]]
    ci = cbind(pred - 2*exp(theta[(length(theta) - ncol(p.sites) + 1):length(theta)][idx.k]),pred + 2*exp(theta[(length(theta) - ncol(p.sites) + 1):length(theta)][idx.k]) )
    plot(dates,obs,main = SiteID[i],
         type="l", ylim = range(ci, na.rm=T))
    polygon(x=c(rev(dates), (dates)), y = c(rev(ci[,1]), (ci[,2])),
            border=NA, col = "grey", density = 100)
    points(dates,pred, type="l", col="red")
    points(y$date[y$ID %in% SiteID[i]],y$obs[y$ID %in% SiteID[i]])
  }

}