R/plot.violinPDFS.R

In mchevalier2/crestr: A Probabilistic Approach to Reconstruct Past Climates Using Palaeoecological Datasets

#' Plot the pdfs as violins
#'
#' Plot the pdfs as violins
#'
#' @inheritParams plot.crestObj
#' @param x A \code{\link{crestObj}} generated by either the \code{\link{crest.calibrate}},
#'        \code{\link{crest.reconstruct}} or \code{\link{crest}} functions.
#' @param climate Climate variables to be used to generate the plot. By default
#'        all the variables are included.
#' @param taxanames A list of taxa to use for the plot (default is all the
#'        recorded taxa).
#' @param ylim The climate range to plot the pdfs on. Default is the full range
#'        used to fit the \code{pdfs} (x$modelling$xrange)
#' @param filename An absolute or relative path that indicates where the diagram
#'        should be saved. Also used to specify the name of the file. Default:
#'        the file is saved in the working directory under the name
#'        \code{'violinPDFs.pdf'}.
#' @param width The width of the output file in inches (default 7.48in ~ 19cm).
#' @param height The height of the output file in inches (default 3in ~ 7.6cm
#'        per variables).
#' @param col A vector of colours that will be linearly interpolated to give a
#'        unique colour to each taxon.
#' @return A table with the climate tolerances of all the taxa
#' @export
#' @examples
#' \dontrun{
#'   data(crest_ex_pse)
#'   data(crest_ex_selection)
#'   reconstr <- crest.get_modern_data(
#'     pse = crest_ex_pse, taxaType = 0,
#'     climate = c("bio1", "bio12"),
#'     selectedTaxa = crest_ex_selection, dbname = "crest_example"
#'   )
#'   reconstr <- crest.calibrate(reconstr,
#'     geoWeighting = TRUE, climateSpaceWeighting = TRUE,
#'     bin_width = c(2, 20), shape = c("normal", "lognormal")
#'   )
#' }
#' ## example using pre-saved reconstruction obtained with the previous command.
#' data(reconstr)
#' ranges <- plot_violinPDFs(reconstr, save=FALSE, ylim=c(5,35),
#'                 taxanames=c(reconstr$inputs$taxa.name[c(2,4,5,1)]),
#'                 col=c('darkblue', 'firebrick3'))
#' lapply(ranges, head)
#'
plot_violinPDFs <- function( x,
                      climate = x$parameters$climate[1],
                      taxanames = x$input$taxa.name,
                      col = viridis::viridis(20), ylim=range(x$modelling$xrange[[climate]]),
                      save = FALSE, filename = 'violinPDFs.pdf',
                      width =  7.48, height = 5,
                      as.png = FALSE, png.res=300
                      ) {

    if(base::missing(x)) x

    if (is.crestObj(x)) {
        test <- is.na(x$modelling$pdfs)
        if( test[1] & length(test) == 1 ) {
            stop('The crestObj requires the climate space to be calibrated. Run crest.calibrate() on your data.\n')
        }

        err <- c()
        for(clim in climate) {
            if(! clim %in% x$parameters$climate) err <- c(err, clim)
        }
        if(length(err) > 0) {
            stop(paste0("The following variables are not available in your crestObj: '", paste(err, collapse="', '"), "'\n\n"))
            return(invisible(NA))
        }

        err <- c()
        for(tax in taxanames) {
            if(! tax %in% x$input$taxa.name) err <- c(err, tax)
        }
        if(length(err) > 0) {
            stop(paste0("The following taxa names are not available in your dataset: '", paste(err, collapse="', '"), "'\n\n"))
            return(invisible(NA))
        }

        if(save) {
            if(as.png) {
                grDevices::png(paste0(strsplit(filename, '.png')[[1]], '.png'), width = width, height = height, units='in', res=png.res)
            } else {
                grDevices::pdf(paste0(strsplit(filename, '.pdf')[[1]], '.pdf'), width=width, height=height)
            }
        } else {
            par_usr <- graphics::par(no.readonly = TRUE)
            on.exit(graphics::par(par_usr))
        }

        ntaxa   <- length(taxanames)
        optima  <- unlist(lapply(x$modelling$pdfs, function(x) return(which.max(x[[climate]]$pdfpol))))
        max_pdf <- unlist(lapply(x$modelling$pdfs, function(x) return(max(x[[climate]]$pdfpol, na.rm=TRUE))))

        if(sum(names(optima) %in% taxanames) ==0) {
            stop("No PDFs were associated with the provided taxa names.\n\n")
        }

        optima <- order(optima[names(optima) %in% taxanames])
        max_pdf <- max_pdf[names(max_pdf) %in% taxanames]


        if(length(taxanames) != length(max_pdf)) {
            warning(paste0("PDFs are not available for all taxa names ('",paste(taxanames[!(taxanames %in% names(max_pdf))], collapse="', '"),"').\n"))
        }
        taxanames <- names(max_pdf)
        ntaxa   <- length(taxanames)

        max_abs <- max(max_pdf, na.rm=TRUE)
        xstep <- mean(max_pdf)/10
        max_pdf <- cumsum(rep(max_pdf[optima] + xstep/2, each=2) / max_abs)
        ystep <- xstep * (diff(ylim) / height) / (max(max_pdf) / width) ## xstep * xyratio

        graphics::par(mar=c(0.5, 2.5, 0.5, 2.5))
        graphics::par(ps=8)

        col <- grDevices::colorRampPalette(col)( ntaxa )

        graphics::par(lwd=0.8)

        graphics::plot(0, 0, type='n',
                       xlim=c(0, max(max_pdf)), ylim=ylim,
                       frame=TRUE, axes=FALSE,
                       xlab='', ylab='', main='', xaxs='i', yaxs='i')

        for(i in graphics::axTicks(2)) {
            if(! i %in% ylim)  graphics::segments(0, i, max(max_pdf), i, col='grey70', lty=3, lwd=0.8)
        }

        for(i in 1:ntaxa) {
            graphics::segments(max_pdf[2*i-1],  ylim[1], max_pdf[2*i-1], ylim[2], col='grey80', lwd=0.5)
            graphics::polygon(c(max_pdf[2*i-1] - x$modelling$pdfs[[taxanames[optima[i]]]][[climate]]$pdfpol / max_abs, rev(max_pdf[2*i-1] + x$modelling$pdfs[[taxanames[optima[i]]]][[climate]]$pdfpol / max_abs)),
                              c(x$modelling$xrange[[climate]], rev(x$modelling$xrange[[climate]])), col=col[i], border='grey80', lwd=0.5)
            graphics::text(max_pdf[2*i-1] - xstep, ylim[1] + ystep, taxanames[optima[i]], srt=90, adj=c(0,0), font=2)
        }

        graphics::par(las=0)
        if(x$misc$dbname == 'private-database') {
            graphics::mtext(climate,side=2,line=1.5, cex=1, font=2)
            graphics::mtext(climate,side=4,line=1.3, cex=1, font=2)
        } else {
            climate_names <- accClimateVariables()
            graphics::mtext(climate_names[climate_names[, 2] == climate, 3],side=2,line=1.5, cex=1, font=2)
            graphics::mtext(climate_names[climate_names[, 2] == climate, 3],side=4,line=1.3, cex=1, font=2)
        }

        graphics::par(mgp=c(2,0.5,0), las=1)
        graphics::axis(2, lwd.ticks=0.8, lwd=0, tck=-0.01, cex.axis=6/7)
        graphics::par(mgp=c(2,0.5,0), las=1)
        graphics::axis(4, lwd.ticks=0.8, lwd=0, tck=-0.01, cex.axis=6/7)

        graphics::rect(0, ylim[1], max(max_pdf), ylim[2], lwd=0.8)

        if(save) {
            grDevices::dev.off()
        }

    } else {
        cat('This function only works with a crestObj.\n\n')
    }
    invisible(pdf_ranges(x, climate=climate, taxanames=taxanames))
}