R/autoplot_ts_stations.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 a timeseries plot with different options for indices and aggregation described below. \cr
#'To see some example outputs of this function, look at the vignette "autoplot-functions" accessible through the package help main page or by typing \cr
#' \emph{vignette("autoplot-functions",package="ClimIndVis")} into the console.
#'
#'@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 one index for one dataset. Returns a seperate plot for each station and each aggregation.
#'   \item "multi_ind" plots several indices for one dataset. Returns a seperate plot for each station and each aggregation. (Ensemble not implemented).
#'   \item "multi_agg" plots 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 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 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.
#'   \item "spi_barplot" plots the SPI value calculated for one station as a barplot. Returns a seperate plot for each station. Different to the other types, the spi_barplots plots the whole time series and does not produce aggregations. For this plot trend calculation is not yet implemented. For a SPI time series of a single month use ts_type="single_ts". For index_args of SPI see \code{\link{index_arguments.spi}}
#'}
#'@param selyears Array to select subset of years to plot, e.g. c(1981:1990).
#'@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.
#'@param ... Optional parameters for shading. Shading only works for type = "single_ts"
#' \itemize{
#' \item \strong{pcols} Array of colors with required length, e.g. 1 for a single timeseries and number of indices, aggregations, points, datasets for <ts_type>. If no colors are provided, colors are set to default colors.
#' \item \strong{shading} Optional. Logical. Add shading of areas for different areas of y-values to plot ? If TRUE nshading ( and yval_shading) need to be provided. See Details for more information.
#'\item \strong{nshading} Optional. Number of shading categories.
#'\item \strong{yval_shading} Optional. Array of the same length as \emph{nshading}+1 with user defined y-limits for shading polygons.
#'\item \strong{col_shading} Optional. Color for shading. Array of same length as nshading.
#'
#' }
#'
#'
#'@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.
#'
#' If you set trendplots = TRUE, a trend line with confidence interval is added to the plot. For every index a default method is used for the trend calculation. The default trend estimation and test depends on the nature of the calculated index:
#' For continuous data (e.g. sum, mean, prcptot, spi) by default an ordinary linear regression is calculated, for count data (e.g dd, cdd, cwd) a logistic regression with time as predictor for the trend estimation and applying a students-t test for significance testing.
#' If you set trend = "MannKendall" inside the list of index arguments (index_args), the non-parametric Mann-Kendall test based on the relative ranking in the series is applied and a Theil-Sen slope is estimated. See e.g. Yue et al. 2002 or Mann (1945), Kendall (1975) for further references.
#'
#'
#'@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_ts_stations(
#'       object_st, index = "sdii",
#'       index_args = list(aggt = "other", aggmons = c(1:3)), trendplots=TRUE,
#'       ts_type = "single_ts", pcols = "royalblue", ylims = c(0,30),
#'       output = "png", plotdir = plotdir, plotname = Sys.Date())
#'
#'  # add shading
#' autoplot_ts_stations(
#'       object_st, index = "tnn",
#'       index_args = list(aggt = "other", aggmons = c(4)), trendplots = TRUE,
#'       ts_type = "single_ts", pcols = "royalblue", ylims = c(-10,20), 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_ts_stations(
#'       dat_p = object_st, index = index, index_args = index_args,
#'       ts_type = "multi_ind", trendplots = TRUE,
#'       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_ts_stations(
#'       object_st, index = list("txn"), index_args=index_args,
#'       ts_type = "multi_agg",
#'       output = "png", plotdir = plotdir)
#'
#'
#' ## ts_typ = "multi_point"
#' autoplot_ts_stations(
#'       object_st, c("dd"), index_args=list(aggt = "other", aggmons = c(1:2)),
#'       ts_type = "multi_point", trendplots = TRUE, selpoints=c(2,3),
#'       output = "png", plotdir = plotdir, plot_title = FALSE, ylims = c(0,100))
#'
#' ## ts_type = "multi_dat"
#'autoplot_ts_stations(
#'      list(object_st, object_hc_st), index = "txn", index_args=list(aggt="monthly", selagg=5),
#'      ts_type = "multi_dat", trendplots = FALSE,
#'      output = "png", plotdir = plotdir, plotname = Sys.Date())
#'
#' ## ts_type = "spi_barplot"
#'autoplot_ts_stations(
#'      object_st, index = "spi", index_args=list(aggt="monthly"),
#'      ts_type = "spi_barplot", plot_title = TRUE)

#'@family autoplot_functions
#'@export
autoplot_ts_stations<-function(dat_p, index, index_args, ts_type = "single_ts",
                               trendplots = FALSE, selpoints=NULL, selyears=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(...)

  # 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" & !all(grepl("spi",index))) stop ("ts_type=spi_barplot only works for index=spi.")
    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)

  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,selyears=selyears)

  # 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()}  )

  # 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_dat, function(x) x$index_info$aggt)))> 1,
    tdims<-get_matching_tdims_index_multi_agg(ind_dat,...),
    tdims<-get_matching_tdims_index(ind_dat,...))
  pdims <- match_pdims_index(ind_dat,selpoints=selpoints,...)
  pinfo <- get_pinfo(ind_dat,ts_type,ylims,opargs)
  
  opargs$units = get_leg_units(ind_dat,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,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))
}