R/mgMap.R

In kkeenan02/memgeneMod: Spatial pattern detection in genetic distance data using Moran's Eigenvector Maps

mgMap <-
function(coords, memgene, wid=NULL, hei=NULL, dev.open=FALSE, add.plot=FALSE, legend=FALSE, ...) {
    
    ## FUNCTIONS FROM NUMERICAL ECOLOGY WITH R (Borcard et al, 2011)
    ## included here to support MEM analyses
    sr.value <- function (dfxy, z, xax = 1, yax = 2, method = c("bubble",
        "greylevel"), zmax = NULL, csize = 1, cpoint = 0, pch = 20,
        clegend = 0.75, neig = NULL, cneig = 1, xlim = NULL, ylim = NULL,
        grid = TRUE, addaxes = TRUE, cgrid = 0.75, include.origin = TRUE,
        origin = c(0, 0), sub = "", csub = 1, possub = "topleft",
        pixmap = NULL, contour = NULL, area = NULL, add.plot = FALSE) {
    
    #
    # Slightly modified version of ade4's s.value() graphical function.
    # Draws round instead of square bubbles in some plots when argument 
    # "bubble" is called.
    #
    # License: GPL-2
    # Author of the original function s.value: Daniel Chessel
    # Modification: Francois Gillet, 25 August 2012
    #
        dfxy <- data.frame(dfxy)
        if (length(z) != nrow(dfxy))
            stop(paste("Non equal row numbers", nrow(dfxy), length(z)))
        opar <- par(mar = par("mar"))
        on.exit(par(opar))
        par(mar = c(0.1, 0.1, 0.1, 0.1))
        coo <- scatterutil.base(dfxy = dfxy, xax = xax, yax = yax,
            xlim = xlim, ylim = ylim, grid = grid, addaxes = addaxes,
            cgrid = cgrid, include.origin = include.origin, origin = origin,
            sub = sub, csub = csub, possub = possub, pixmap = pixmap,
            contour = contour, area = area, add.plot = add.plot)
        if (!is.null(neig))
        {
            if (is.null(class(neig))) neig <- NULL
            if (class(neig) != "neig") neig <- NULL
            deg <- attr(neig, "degrees")
            if (length(deg) != length(coo$x)) neig <- NULL
        }
        if (!is.null(neig))
        {
            fun <- function(x, coo)
            {
                segments(coo$x[x[1]], coo$y[x[1]], coo$x[x[2]], coo$y[x[2]],
                    lwd = par("lwd") * cneig)
            }
            apply(unclass(neig), 1, fun, coo = coo)
        }
        method <- method[1]
        if (method == "greylevel")
        {
            br0 <- pretty(z, 6)
            nborn <- length(br0)
            coeff <- diff(par("usr")[1:2])/15
            numclass <- cut.default(z, br0, include.lowest = TRUE, labels = FALSE)
            valgris <- seq(1, 0, le = (nborn - 1))
            h <- csize * coeff
            for (i in 1:(nrow(dfxy)))
                {
                    symbols(coo$x[i], coo$y[i], circles = h/2, 
                        bg = gray(valgris[numclass[i]]),
                        add = TRUE, inches = FALSE)
                }
            scatterutil.legend.circle.grey(br0, valgris, h/2, clegend)
            if (cpoint > 0) points(coo$x, coo$y, pch = pch, cex = par("cex") * cpoint)
        }
        else if (method == "bubble")
        {
            coeff <- diff(par("usr")[1:2])/15
            sq <- sqrt(abs(z))
            if (is.null(zmax)) zmax <- max(abs(z))
            w1 <- sqrt(zmax)
            sq <- csize * coeff * sq/w1
            for (i in 1:(nrow(dfxy)))
            {
                if (sign(z[i]) >= 0)
                {
                    symbols(coo$x[i], coo$y[i], circles = sq[i]/2, bg = "black", 
                        fg = "white", add = TRUE, inches = FALSE)
                }
                else
                {
                    symbols(coo$x[i], coo$y[i], circles = sq[i]/2, bg = "white", 
                        fg = "black", add = TRUE, inches = FALSE)
                }
            }
            br0 <- pretty(z, 4)
            l0 <- length(br0)
            br0 <- (br0[1:(l0 - 1)] + br0[2:l0])/2
            sq0 <- sqrt(abs(br0))
            sq0 <- csize * coeff * sq0/w1
            sig0 <- sign(br0)
            if (clegend > 0) scatterutil.legend.bw.circle(br0, sq0, sig0, clegend)
            if (cpoint > 0) points(coo$x, coo$y, pch = pch, cex = par("cex") * cpoint)
        }
        else if (method == "circlesize") print("not yet implemented")
        if (!add.plot) box()
        invisible(match.call())
    }
    
    
    
    scatterutil.legend.bw.circle <- function (br0, sq0, sig0, clegend) {
        br0 <- round(br0, digits = 6)
        cha <- as.character(br0[1])
        for (i in (2:(length(br0)))) cha <- paste(cha, br0[i], sep = " ")
        cex0 <- par("cex") * clegend
        yh <- max(c(strheight(cha, cex = cex0), sq0))
        h <- strheight(cha, cex = cex0)
        y0 <- par("usr")[3] + yh/2 + h/2
        ltot <- strwidth(cha, cex = cex0) + sum(sq0) + h
        rect(par("usr")[1] + h/4, y0 - yh/2 - h/4, 
            par("usr")[1] + ltot + h/4, y0 + yh/2 + h/4, col = "white")
        x0 <- par("usr")[1] + h/2
        for (i in (1:(length(sq0))))
        {
            cha <- br0[i]
            cha <- paste(" ", cha, sep = "")
            xh <- strwidth(cha, cex = cex0)
            text(x0 + xh/2, y0, cha, cex = cex0)
            z0 <- sq0[i]
            x0 <- x0 + xh + z0/2
            if (sig0[i] >= 0)
                symbols(x0, y0, circles = z0/2, bg = "black", fg = "white",
                    add = TRUE, inches = FALSE)
            else symbols(x0, y0, circles = z0/2, bg = "white", fg = "black",
                add = TRUE, inches = FALSE)
            x0 <- x0 + z0/2
        }
        invisible()
    }
    
    
    
    scatterutil.legend.circle.grey <- function (br0, valgris, h, clegend) {
        if (clegend <= 0) return(invisible())
        br0 <- round(br0, digits = 6)
        nborn <- length(br0)
        cex0 <- par("cex") * clegend
        x0 <- par("usr")[1] + h
        x1 <- x0
        for (i in (2:(nborn)))
        {
            x1 <- x1 + h
            cha <- br0[i]
            cha <- paste(cha, "]", sep = "")
            xh <- strwidth(cha, cex = cex0)
            if (i == (nborn)) break
            x1 <- x1 + xh + h
        }
        yh <- max(strheight(paste(br0), cex = cex0), h)
        y0 <- par("usr")[3] + yh/2 + h/2
        rect(par("usr")[1] + h/4, y0 - yh/2 - h/4, x1 - h/4, y0 + yh/2 + h/4, 
            col = "white")
        x0 <- par("usr")[1] + h
        for (i in (2:(nborn)))
        {
            symbols(x0, y0, circles = h/2, bg = gray(valgris[i - 1]), add = TRUE, 
                inches = FALSE)
            x0 <- x0 + h
            cha <- br0[i]
            if (cha < 1e-05) cha <- round(cha, digits = 3)
            cha <- paste(cha, "]", sep = "")
            xh <- strwidth(cha, cex = cex0)
            if (i == (nborn)) break
            text(x0 + xh/2, y0, cha, cex = cex0)
            x0 <- x0 + xh + h
        }
        invisible()
    }
    
    
    
    
    
    ## mgMap specific code
    
    memgene <- as.matrix(memgene)
    
    if (!legend) {
        clegend <- -1
    }
    else {
        clegend <- 0.75
    }
    
    if (is.null(wid) && is.null(hei) && (ncol(memgene)==2) && !add.plot) {
        sideBySide <- TRUE
    }
    else {
        sideBySide <- FALSE
    }

    if (sideBySide) {
        if (!dev.open) {
            dev.new(width=9, height=4.5)
        }
        par(mfcol=c(1,2))
        par(mar=c(2,2,2,2))
        plot(coords, type="n", main="", xlab="", ylab="")
        sr.value(coords, memgene[,1], add.plot=TRUE, clegend=clegend, ...)
        plot(coords, type="n", main="", xlab="", ylab="")
        sr.value(coords, memgene[,2], add.plot=TRUE, clegend=clegend, ...)
    }
    else {
        if (is.null(wid)) {
            wid <- 7
            hei <- 7
        }
        for (i in 1:ncol(memgene)) {
            if (!add.plot) {
                if (!dev.open) {
                    dev.new(width=wid, height=hei)
                }
                par(mar=c(2,2,2,2))
                plot(coords, type="n", main="", xlab="", ylab="")
                sr.value(coords, memgene[,i], add.plot=TRUE, clegend=clegend, ...)

            }
            else {
                sr.value(coords, memgene[,i], add.plot=TRUE, clegend=clegend, ...)
            }
        }
    }
}