R/plot.scatterPDFS.R

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

#' Plot the pdf optima and uncertainty ranges in a climate biplot
#'
#' Plot the pdf optima and uncertainty ranges in a climate biplot
#'
#' @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 Names of the two climate variables to be used to generate the plot. By default
#'        plot. By default the first two variables are included.
#' @param taxanames A list of taxa to use for the plot (default is all the
#'        recorded taxa).
#' @param xlim, ylim The climate range to plot the data. Default is the full range
#'        of the observed climate space.
#' @param add_modern A boolean to add the location and the modern climate values
#'        to the plot (default \code{FALSE}).
#' @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).
#' @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)
#' dat <- plot_scatterPDFs(reconstr, save=FALSE,
#'                  taxanames=c(reconstr$inputs$taxa.name[c(2,4,5,1)]))
#' dat
#'
plot_scatterPDFs <- function( x,
                      climate = x$parameters$climate[1:2],
                      taxanames = x$input$taxa.name,
                      uncertainties = x$parameters$uncertainties,
                      xlim=range(x$modelling$climate_space[, climate[1]]),
                      ylim=range(x$modelling$climate_space[, climate[2]]),
                      add_modern = FALSE,
                      save = FALSE, filename = 'scatterPDFs.pdf',
                      width = 5.51, height = 5.51,
                      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))
        }

        COL1 = c('#d95f02', '#1b9e77', '#7570b3', 'white')

        get_col <- function(use) {
            use <- x$inputs$selectedTaxa[use, climate[1]] + 10*x$inputs$selectedTaxa[use, climate[2]]
            if(use == 1) return(COL1[1])
            if(use == 10) return(COL1[2])
            if(use == 11) return(COL1[3])
            return(COL1[4])
        }

        if(add_modern) {
            if (length(x$misc$site_info) <= 3) {
                add_modern <- FALSE
            }
        }

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

        climate <- climate[climate %in% x$parameters$climate]
        if(length(climate) > 2) {
            climate <- climate[1:2]
            warning('Only the first two valid variable names were used for the plot.\n\n')
        }
        if(length(climate) < 2) stop('Two valid climate variable names are necessary.\n\n')

        ntaxa   <- length(taxanames)
        optima1  <- unlist(lapply(x$modelling$pdfs, function(x) return(which.max(x[[climate[1]]]$pdfpol))))
        optima2  <- unlist(lapply(x$modelling$pdfs, function(x) return(which.max(x[[climate[2]]]$pdfpol))))

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

        n <- names(optima1)
        optima1 <- x$modelling$xrange[[climate[1]]][optima1]  ;  names(optima1) <- n
        optima2 <- x$modelling$xrange[[climate[2]]][optima2]  ;  names(optima2) <- n

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

        ntaxa   <- length(taxanames)
        ranges <- pdf_ranges(x, climate=climate, taxanames=taxanames, uncertainties=sort(uncertainties, decreasing=TRUE))

        graphics::layout(matrix(c(2, 1), ncol=1), height=c(0.08, 0.92))

        {
            graphics::par(mar=c(2.5, 2.5, 0, 0.5))
            graphics::par(ps=8)
            graphics::par(lwd=0.8)


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

            graphics::abline(v = graphics::axTicks(1), col='grey70')
            graphics::abline(h = graphics::axTicks(2), col='grey70')
            graphics::points(x$modelling$climate_space[, climate[1]], x$modelling$climate_space[, climate[2]], pch=15, cex=0.5, col='grey40')

            if(add_modern) {
                if(is.numeric(x$misc$site_info$climate[, climate[1]])) {
                    graphics::points(x$misc$site_info$climate[, climate[1]], x$misc$site_info$climate[, climate[2]], pch=24, col=NA, bg='red', cex=1.5, lwd=1.5)
                }
            }

            for(tax in taxanames) {
                for(u in 1:length(uncertainties)) {
                    graphics::segments(ranges[[u]][tax, 1], optima2[tax], ranges[[u]][tax, 2], optima2[tax], lwd=u/length(uncertainties))
                    graphics::segments(optima1[[tax]], ranges[[u]][tax, 4], optima1[[tax]], ranges[[u]][tax, 5], lwd=u/length(uncertainties))
                }
                graphics::points(optima1[tax], optima2[tax], pch=21, cex=1.5, bg=sapply(names(optima1[tax]), get_col, simplify=TRUE), col='black', lwd=1.2)

            }

            graphics::par(mgp=c(2,0.3,0), las=1)
            graphics::axis(2, lwd.ticks=0, lwd=0, tck=-0, cex.axis=7/8)
            graphics::par(mgp=c(2,0,0), las=1)
            graphics::axis(1, lwd.ticks=0, lwd=0, tck=0, cex.axis=7/8)


            climate_names <- accClimateVariables()
            graphics::mtext(paste0(climate_names[climate_names[, 2] == climate[1], 3], ' [', climate[1], ']'), side=1, line=1, cex=1, font=2)
            graphics::par(las=0)
            graphics::mtext(paste0(climate_names[climate_names[, 2] == climate[2], 3], ' [', climate[2], ']'), side=2, line=1.3, cex=1, font=2)
            graphics::par(las=1)
        }

        {
            graphics::par(mar=c(0, 2.5, 0, 0.5))
            graphics::plot(0, 0, type='n', xlim=c(0,1), ylim=c(0, 1), xaxs='i', yaxs='i', frame=FALSE, axes=FALSE, xlab='', ylab='')
            graphics::points(c(0.02, 0.27, 0.52, 0.77), rep(0.5, 4), col='black', bg=COL1, pch=21, cex=2, lwd=1.2)
            graphics::text(0.05, 0.5, paste0(climate[1], ' only'), cex=8/8, adj=c(0, 0.5))
            graphics::text(0.30, 0.5, paste0(climate[2], ' only'), cex=8/8, adj=c(0, 0.5))
            graphics::text(0.55, 0.5, 'Both variables', cex=9/8, adj=c(0, 0.5))
            graphics::text(0.80, 0.5, 'Not selected', cex=9/8, adj=c(0, 0.5))
        }

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

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