R/autoplot_anomaly_ts.R

In Climandes/ClimIndVis: Calculation, verification and visualisation of climate indices

#'Wrapper function to calculate and plot timerseries of indices from station data
#'
#'The function calls \code{\link{calc_index}} and returns an anomaly timeseries plot with different options for indices and aggregation described below. \cr
#'At the moment the function does not have the options to plot trends, please use autoplot_ts_stations for plotting trends. \cr
#'
#'
#'@inheritParams plot_doc
#'@param dat_p A Climindvis-object  with type = "p" generated by \code{\link{make_object}}. For ts_type="multi_dat" dat_p is a list of climindivs-objects with type = "p" or "p_hc".
#'@param index_args List of arguments for index. See \code{\link{calc_index}}.\cr For ts_type = "multi_ind", index_args is a list containing a list of index arguments for each of the indices.\cr For ts_type="multi_agg" index_args is a list containing a list of index arguments for each aggregation type (e.g. if you want to compare annual and seasonal values type index_args=list(list(aggt="annual"),list(aggt="seasonal")).
#'@param index Name of index to be calculated or array of indices to be calculated for ts_type = "multi_ind".
#'@param ts_type Type of time series plot to be drawn:
#' \itemize{
#'   \item "single_ts" plots the anomaly of one index for one dataset. Returns a seperate plot for each station and each aggregation.
#'   \item "multi_ind" plots the anomalies of several indices for one dataset. Returns a seperate plot for each station and each aggregation. (Ensemble not implemented).
#'   \item "multi_agg" plots the anomaly of one index and one dataset. Returns a seperate plot for each station and combines aggregation values in one plot. (Ensemble not implemented).
#'   \item "multi_point" plots the anomaly of one index and one dataset. Returns a seperate plot for each aggregation value and combines several stations in one plot (either all stations or a subset to be specified in selpoints). (Ensemble not implemented).
#'   \item "multi_dat" plots the anomaly of one index for multiple datasets (e.g. observations, reanalyses and hindcasts). Returns a seperate plot for each station and each aggregation. If ensemble data is provided, the spread of the ensemble members is plotted.
#'}
#'@param refyears Integer array of reference years for calculation of anomaly. Default=1981:2010
#'@param ylims Limits of xaxis for plot of format ylims=c(ymin,ymax). If a individual limits for every plot are needed, set ylims to "ind". If not provided limits are taken from the data and are equal in all plots.
#'@param plot_title Logical. Add plot title and additional information to graphic? Default = TRUE.
#'@param plot_legend Logical. Add plot legend of depicted stations to graphic? Default = TRUE.
#'@param cex cex of plots (see par()). Default=1
#'@param text_cex size of text plotted in graph. Default=1
#'@param pwidth Width of plotting region of output (as input for jpeg, png, ...). For output = pdf, pwidth is set to pwidth/pres. Default = 2000.
#'@param pheight Height of plotting region of output (as input for jpeg, png, ...). For output = pdf, pheight is set to pheight/pres.  Default = 1300.
#'@param pres Graphics resolution. Default = 250.
#'
#'
#'
#'@section Details:
#'
#' The missing values displayed when plot_title = TRUE denotes the percentage of missing daily values in the selected aggregation period.
#'
#' By passing the optional parameters \emph{shading}, \emph{nshading} (\emph{ycal_shading}) and \emph{col_shading} to the function, horizontal shading polygons can be added to the plot. \cr
#' For example, to mark the area below zero in blue and the one above zero in red when looking at temperature indices such as maximum teperatures one has to pass the following values to the function: \emph{shading}=TRUE, \emph{nshading}=2, \emph{yval_shading}=c(-999,0,999), \emph{col_shading}=c("blue","red").\cr
#' If you want to add equidistant shading categories (e.g. if you issue a tercile forecast and want to show the corresponding values in the climatology you just need to specify \emph{nshading} which would be equal to the number of forecast categories and the function calculates the corresponding category thresholds. This is used e.g. by the \code{\link{autoplot_forecast_stations}} function when \emph{plot_climatology}=TRUE.
#'
#' Trendplots are currently not implemented in this function.
#'
#'@examples
#'  data("object_st","object_hc_st") # load example data
#' plotdir="./" # define a plot directory here or change output to NULL or "dev.new".
#'
#' ## ts_type = single_ts
#' autoplot_anomaly_ts(
#'       object_st, index = "sdii",
#'       index_args = list(aggt = "other", aggmons = c(1:3)),
#'       ts_type = "single_ts", pcols = "royalblue",
#'       output = "png", plotdir = plotdir, plotname = Sys.Date())
#'
#'  # add shading
#' autoplot_anomaly_ts(
#'       object_st, index = "tnn",
#'       index_args = list(aggt = "other", aggmons = c(4)),
#'       ts_type = "single_ts", pcols = "royalblue", selpoints=2,
#'       shading = TRUE, nshading = 2, yval_shading = c(-999,0,999),col_shading = c("blue","red"),
#'       output = "png", plotdir = plotdir, plotname = Sys.Date())
#'
#' ## ts_type = multi_ind
#' index = c("txn", "txx","dd")
#' index_args = list(
#'       txn = list(aggt = "other", aggmons = c(1:3)),
#'       txx = list(aggt = "other", aggmons = c(1:3)),
#'       tnx = list(aggt = "other", aggmons = c(1:3)))
#' autoplot_anomaly_ts(
#'       dat_p = object_st, index = index, index_args = index_args,
#'       ts_type = "multi_ind",
#'       output = "png", plotdir = plotdir, plotname = Sys.Date())
#'
#'
#' ## ts_type = "multi_agg"
#' index_args = list (
#'       list(aggt="seasonal",selagg=c("MAM")),
#'       list(aggt="monthly", selagg=c("1","2","3")))
#' autoplot_anomaly_ts(
#'       object_st, index = list("txn"), index_args=index_args,
#'       ts_type = "multi_agg",
#'       output = "png", plotdir = plotdir)
#'
#'
#' ## ts_typ = "multi_point"
#' autoplot_anomaly_ts(
#'       object_st, c("dd"), index_args=list(aggt = "other", aggmons = c(1:2)),
#'       ts_type = "multi_point", selpoints=c(2,3),
#'       output = "png", plotdir = plotdir, plot_title = FALSE)
#'
#' ## ts_type = "multi_dat"
#'autoplot_anomaly_ts(
#'      list(object_st, object_hc_st), index = "txn", index_args=list(aggt="monthly", selagg=5),
#'      ts_type = "multi_dat",
#'      output = "png", plotdir = plotdir, plotname = Sys.Date())
#'@family autoplot_functions
#'@export
autoplot_anomaly_ts<-function(dat_p, index, index_args, ts_type = "single_ts",
                               selpoints=NULL, refyears = NULL,
                              output = NULL, plotdir = NULL, plotname = NULL, title = "", plot_title = TRUE,plot_legend=TRUE,
                              ylims = NULL, pwidth = 2000, pheight = 1300, pres = 250, cex = 1, text_cex = 1, ...){
  opargs <- list(...)
  trendplots=FALSE # trends not implemented yet
  # 1. checks general --------------------------------------

  points <- ifelse (missing(dat_p), 0, 1)
  grid <- 0

  #multi_dat checks:
  tryCatch({
    if (!any(is.element(ts_type,c("single_ts","multi_ind","multi_agg","multi_point","multi_dat","spi_barplot")))) stop("<<ts_type>> passed to function is not defined. ")
    if (ts_type=="spi_barplot") stop ("ts_type=spi_barplot can only be plottet with autoplot_ts_stations.")
    if (is.null(dat_p$data) && length(dat_p)>1 && !ts_type=="multi_dat"){stop("Choose <multi_dat> for several datasets.")}
  if(ts_type=="multi_dat") {
    if (class(dat_p[[1]])!="climindvis" | length(dat_p) ==1) stop("for ts_type=multi_dat you need to enter a list of climindvis objects")
    dat <- dat_p
    } else {
    dat=get_data_list(c("dat"))
  }}, error=function(cond){
    message("Error in input arguments:")
    message(cond)
    stop_quietly()}  )

  check_dir(plotdir,output)


  #if (length(index)==1  &&  length(index_args)!=1 && !is.list(index_args[[1]])){index_args <- list(index_args)}
  check_class(dat[!sapply(dat,is.null)],"climindvis")
  check_type(dat, index, index_args, ts_type)
  check_dat(dat)

  # trends currently not implemented
  if (is.list(index_args[[1]])){
  index_args <-  lapply(index_args, function(ii) {
    if (ii$aggt=="xdays") stop("aggt=xdays is not implemented yet for autoplot functions")
     ii$trend <- trendplots
     return(ii)})
  } else index_args$trend  <- trendplots

  dat<-cut_to_same_dates(dat)
  datyears=as.integer(levels(dat[[1]]$time_factors$years))
  if (is.null(refyears)) {
    refyears=1981:2010
  } else {
    refyears=as.integer(refyears)
    if( !all(is.element(refyears,datyears))) stop("at least one year in refyears is missing in input data")
  }
  refy=switch(all(diff(refyears)==1)+1,paste(refyears,collapse="_"),paste(range(refyears),collapse="-"))


  # 2. calculations -------------------------------------------
  # 2.1 calculation of index and trend
  tryCatch({ if(length(index) == 1 & is.list(index_args[[1]])) { # multi_agg
    ind_dat<- lapply(1:length(index_args), function(ii)
      do.call("calc_index", c(list(dat[[1]],index[1]),index_args[[ii]])))
  }else if(length(index) == 1){ # single + multi_dat
    ind_dat<-lapply(dat, function(x) do.call("calc_index",c(list(x,index[[1]]),index_args)))

    #ind_dat<-lapply(dat, function(x) calc_index(x,index[[1]],unlist(index_args)))
  } else if(length(dat)== 1 & length(index)>1 ) {  # multi_ind
    ind_dat<- lapply(1:length(index), function(ii)
      do.call("calc_index", c(list(dat[[1]],index[ii]),index_args[[ii]])))
  } else stop("please check your input. Function does not work for multiple climindvis objects and multiple index")
   }, error=function(cond){
    message("Error in input arguments:")
    message(cond)
    stop_quietly()}  )


  ind_anom = lapply(ind_dat, function (ii){
    anom=ii[which(names(ii)!="index")]
  iref=index_array(ii$index,dim=which(ii$index_info$idims=="year"),which(is.element(as.integer(ii$index_info$years),refyears)),drop=FALSE)
  irefm=apply(iref,c(1:length(ii$index_info$idims))[-which(is.element(ii$index_info$idims, c("year", "ens")))],
              mean, na.rm = TRUE)
  # ian=index_array(ii$index,dim=which(ii$index_info$idims=="year"),which(is.element(as.integer(ii$index_info$years),anomyears)),drop=FALSE)
  # ianm=apply(ian,c(1:length(ii$index_info$idims))[-which(is.element(ii$index_info$idims, c("year", "ens")))],
  #            mean, na.rm = TRUE)
  irefm_help <- rep(irefm, length(ii$index_info$years))
  anom$index = ii$index-irefm_help
  return(anom)})


  # 2.2 calculate missing values
  if(!is.list(index_args[[1]])) {index_args <- list(index_args)}

  index2 <- lapply(index, function(x) {x})
  names(index2) <- index

  ivar <- mapply(function(args,index){
    ivar<-get_indexvar(index,args)
    return(list(ivar))
  },args=index_args,index=index2)


  tfact <- mapply(function(args,index,obj){
    class(obj)  <- append(index,"climindvis")
    arguments=c(list(climindvis=obj),args=args)
    tfact<-get_date_factors(obj,args$aggt,args$aggmons,args$selagg,args$start_days,args$end_days)$tfactor
    return(list(tfact))
  },index=index2,args=index_args,obj=dat)


  day_dat_na <- mapply(function(ivar,tfact,dat) {
    if(dat$data_info$type=="p_hc"){
      lev <- unique(gsub("....-","",tfact))
      ll <- length(lev[-is.na(lev)])
      dd <- dim(dat$data[[ivar]])[1]
      na_year <- matrix(rep(0, ifelse(ll>0,ll,1)*dd),ncol=ll, nrow=dd )
    } else {
      var <- dat$data[[ivar]]
      na_var <- apply(var,c(1), function(x) climdex.pcic:::tapply.fast(x, tfact, function(t) sum(is.na(t))*100/length(t)))
      na_agg <- rearange_by_year(na_var,levels(tfact))
      d <- dim(na_agg)
      na_year <-apply(na_agg,1:(length(d)-1), function(x) {
        sum(x, na.rm = TRUE)*100/(length(x)*100)
      })
    }

    return(list(na_year))
  },ivar=ivar,tfact=tfact,dat=dat)


  # 3. Plot
  # 3.1 get plot arguments
  ifelse(ts_type =="multi_agg" &  length(Reduce(union,lapply(ind_anom, function(x) x$index_info$aggt)))> 1,
    tdims<-get_matching_tdims_index_multi_agg(ind_anom,...),
    tdims<-get_matching_tdims_index(ind_anom,...))
    #xxx in ind_anom[[1]]$index_info gibt es $years und $aggnames. Wenn man oben mit cut_to_same_dates ausschneidet braucht man diese Funktione hier gar nicht mehr, oder?
  pdims <- match_pdims_index(ind_anom,selpoints=selpoints,...)
  pinfo <- get_pinfo(ind_anom,ts_type,ylims,opargs)
  opargs$units = get_leg_units(ind_anom,trend=FALSE) #false because we do not want brackets

  if (is.null(selpoints)) selpoints=1:pdims$pnumber

  # 3.2 call respective plotting function
 do.call(ts_type, list(ind_dat=ind_anom,day_dat_na,tdims, pdims, pinfo=pinfo, trendplots, output=output,plotdir=plotdir,plotname=plotname,
                    pwidth=pwidth,pheight=pheight,pres=pres,selpoints=selpoints,cex=cex,text_cex=text_cex,plot_title=plot_title,title=title,plot_legend=plot_legend,opargs=opargs,anom=TRUE,refyears=refy))
}