generic geoplotting functions

geoplots.R

# load deps
library(rgdal)
library(raster)
library(data.table)
library(rgeos)
library(tools)
library(RColorBrewer)
library(dplyr)
library(shape)

plotRaster=function(x, col, add = FALSE, legend = TRUE, horizontal = FALSE, 
    legend.shrink = 0.5, legend.width = 0.6, legend.mar = ifelse(horizontal, 
        3.1, 5.1), legend.lab = NULL, graphics.reset = FALSE, 
    bigplot = NULL, smallplot = NULL, legend.only = FALSE, lab.breaks = NULL, 
    axis.args = NULL, legend.args = NULL, interpolate = FALSE, 
    box = TRUE, breaks = NULL, zlim = NULL, zlimcol = NULL, fun = NULL, 
    asp, colNA = NA, ...) {

    ffun <- NULL
    if (is.character(fun)) {
        if (fun %in% c("sqrt", "log")) {
            if (fun == "sqrt") {
                ffun <- fun
                fun <- sqrt
            }
            else {
                ffun <- fun
                fun <- log
            }
        }
        else {
            fun - NULL
        }
    }
    else {
        fun <- NULL
    }
    if (missing(asp)) {
        if (isLonLat(x)) {
            ym <- mean(c(x@extent@ymax, x@extent@ymin))
            asp <- 1/cos((ym * pi)/180)
        }
        else {
            asp <- 1
        }
    }
    e <- as.vector(t(bbox(extent(x))))
    x <- as.matrix(x)
    if (!is.null(fun)) {
        x <- fun(x)
    }
    x[is.infinite(x)] <- NA
    if (!is.null(zlim)) {
        if (!is.null(zlimcol)) {
            x[x < zlim[1]] <- zlim[1]
            x[x > zlim[2]] <- zlim[2]
        }
        else {
            x[x < zlim[1] | x > zlim[2]] <- NA.com/
        }
    }
    w <- getOption("warn")
    options(warn = -1)
    if (is.null(breaks)) {
        zrange <- range(x, zlim, na.rm = TRUE)
    }
    else {
        zrange <- range(x, zlim, breaks, na.rm = TRUE)
    }
    options(warn = w)
    if (!is.finite(zrange[1])) {
        legend <- FALSE
    }
    else {
        x <- raster:::.asRaster(x, col, breaks, zrange, colNA)
    }
    old.par <- par(no.readonly = TRUE)
    if (add) {
        big.plot <- old.par$plt
    }
    if (legend.only) {
        graphics.reset <- TRUE
    }
    if (is.null(legend.mar)) {
        legend.mar <- ifelse(horizontal, 3.1, 5.1)
    }
    # temp <- .imageplotplt(add = add, legend.shrink = legend.shrink, 
        # legend.width = legend.width, legend.mar = legend.mar, 
        # horizontal = horizontal, bigplot = bigplot, smallplot = smallplot)
    # smallplot <- temp$smallplot
    # bigplot <- temp$bigplot
	rasterImage(x, e[1], e[3], e[2], e[4], interpolate = interpolate)
}


singleMapFun=function(rast, studyareaSHP=NULL, basemapSHP=NULL, pointsSHP=NULL, linesSHP=NULL, palname=NULL, revpal=FALSE, digit=NULL, cols=NULL, centreCols=FALSE, main=NULL, path=NULL, ncols=100, maxpixels=ncell(rast), res="low", letter=NULL, colScale="quantile") {
	# generate colors
	if (!is.null(palname)) {
		cols=brewer.pal(brewer.pal.info[palname,"maxcolors"], palname)
		if (revpal) {
			cols=rev(cols)
		}
	} else if (!is.null(cols)) {
		cols=rgb(colorRamp(c("blue", "red", "orange", "yellow"))(seq(0, 1, length.out=ncols)), maxColorValue=255)
	} else {
		cols=rgb(colorRamp(c("white", "black"))(seq(0, 1, length.out=ncols)), maxColorValue=255)
	}
	
	# calculate color breaks
	if (colScale=="quantile") {
		# generate quantile-based brks
		brks=as.numeric(quantile(rast, probs=seq(0,1,length.out=length(cols))))
		for (i in seq_along(brks)[-1]) {
			if (brks[i]==brks[1]) {
				brks[i]=brks[i-1]+1e-05
			}
		}
	} else {
		# generate linear-based brks
		if (!centreCols) {
			brks=seq(cellStats(rast, "min"), cellStats(rast, "max"), length.out=length(cols))
		} else {
			mid=ceiling(length(cols)/2)
			first=1
			last=length(cols)
			brks=numeric(length(cols))
			brks[first:mid]=seq(cellStats(rast, "min"), 0, length.out=length(first:mid))
			brks[mid:last]=seq(0, cellStats(rast, "max"), length.out=length(mid:last))		
		}
	}

	# set parameters if plot is saved at high or low resolution
	if (res=="high") {
		par(cex.lab=0.5, cex.main=1)
		res=500
		height=8
		width=8
		units="cm"
	} else {
		par(cex.lab=0.5, cex.main=1)
		res=100
		height=8
		width=8
		units="cm"
	}
	
	# sink plot to file if path supplied
	if (!is.null(path)) {
		do.call(file_ext(path), args=list(file=path, height=height, width=width, res=res, units=units))
	}

	if (!is.null(main) && nchar(main)>0) {
		# set plotting window
		layout(matrix(c(1,1,3,2),ncol=2,nrow=2,byrow=TRUE), height=c(0.1, 0.9), width=c(0.9,0.3))

		# make figure title
		par(mar=c(0.1,0.1,0.1,0.1))
		plot(1, xlab="", ylab="", main="", bty="n", type="n", axes=FALSE)
		mtext(side=1, at=1, text=main, line=-1)
	} else {
		# set plotting window
		layout(matrix(c(2,1),ncol=2,nrow=1,byrow=TRUE), height=c(1, 1), width=c(0.9,0.3))	
	}
	
	# make color legend
	par(mar=c(0.1,0.1,0.1,0.1))
	plot(1, xlab="", ylab="", main="", bty="n", type="n", axes=FALSE)
	if (is.null(letter)) {
		posy=c(0.1, 0.9)
	} else {
		posy=c(0.05, 0.8)
	}
	if(is.null(digit))
		digit=ifelse(maxval>1, 0, 1)
	if (colScale=='linear') {
		legValues=pretty(brks)
		legValues=legValues[which(legValues > min(brks) & legValues < max(brks))]
		shape::colorlegend(zlim=c(min(brks), max(brks)), digit=digit, col=cols, zval=legValues, posx=c(0.3, 0.4), posy = posy, xpd=TRUE)
	} else {
		legValues=pretty(seq(0,1,length.out=length(cols)))
		legValues=legValues[which(legValues >= 0 & legValues <= 1)]
		labValues=as.numeric(quantile(rast, probs=legValues))
		colorlegend(zlim=c(0,1), digit=digit, col=cols, zval=legValues, zlabel=labValues, posx=c(0.3, 0.4), posy = posy , xpd=TRUE)
	}
	# add letter
	if (!is.null(letter))
		mtext(side=3, text=paste0('(', letter, ')'), line=-2, at=par()$usr[1]+(abs(diff(par()$usr[1:2]))*0.325), cex=1.5)
		
	
	# geoplot
	par(mar=c(0.5,0.5,0.5,0.5))
	if (!is.null(studyareaSHP)) {
		xlimit=studyareaSHP@bbox[1,]
		ylimit=studyareaSHP@bbox[2,]
	} else {
		xlimit=c(xmin(rast), xmax(rast))
		ylimit=c(ymin(rast), ymax(rast))
	}
	plot(1, xlim=xlimit, ylim=ylimit, type="n", axes=FALSE, bty="n")
	if (!is.null(basemapSHP)) {
		plot(basemapSHP, col="grey80", border="grey90", add=TRUE)
	}
	plotRaster(rast, col=cols, maxpixels=maxpixels, breaks=brks, legend=FALSE, box=FALSE, axes=FALSE,  bty="n", main=main, add=TRUE)
	if (!is.null(linesSHP)) {
		lines(linesSHP, lwd=1.5, col='black')
	}
	if (!is.null(pointsSHP)) {
		points(pointsSHP, col='black', cex=1.2, pch=16)
	}
	if (!is.null(basemapSHP)) {
		plot(basemapSHP, border="black", add=TRUE)
	}

	# close plotting device if saving to file
	if (!is.null(path)) {
		dev.off()
	}
}



dualMapFun=function(rast1, rast2, studyareaSHP=NULL, basemapSHP=NULL, palname=NULL, revpal=FALSE, digit=NULL, cols=NULL, centreCols=FALSE, main=NULL, path=NULL, res="high", ncols=100, maxpixels=100000, ptsLIST=NULL, ptsCols=NULL, letter=NULL) {
	# generate colors
	if (!is.null(palname)) {
		cols=brewer.pal(brewer.pal.info[palname,"maxcolors"], palname)
		if (revpal) {
			cols=rev(cols)
		}
	} else if (!is.null(cols)) {
		cols=rgb(colorRamp(cols)(seq(0, 1, length.out=ncols)), maxColorValue=255)
	} else {
		cols=rgb(colorRamp(c("white", "black"))(seq(0, 1, length.out=ncols)), maxColorValue=255)
	}
	
	# generate brks using linear color scale
	maxval=max(cellStats(rast1, "max"), cellStats(rast2, "max"))
	minval=max(cellStats(rast1, "min"), cellStats(rast2, "min"))	
	if (!centreCols) {
		brks=seq(minval, maxval, length.out=length(cols))
	} else {
		mid=ceiling(length(cols)/2)
		first=1
		last=length(cols)
		brks=numeric(length(cols))
		brks[first:mid]=seq(minval, 0, length.out=length(first:mid))
		brks[mid:last]=seq(0, maxval, length.out=length(mid:last))
	}

	# set parameters if plot is saved at high or low resolution
	if (res=="high") {
		par(cex.lab=0.5, cex.main=1)
		res=500
		height=8
		width=17
		units="cm"
	} else {
		par(cex.lab=0.5, cex.main=1)
		res=100
		height=8
		width=8
		units="cm"		
	}
	
	# sink plot to file if path supplied
	if (!is.null(path)) {
		do.call(file_ext(path), args=list(file=path, height=height, width=width, res=res, units=units))
	}

	# set plotting window
	layout(matrix(c(2,3,1),ncol=3,nrow=1,byrow=TRUE), width=c(1,1,0.3))
	
	# make color legend
	par(mar=c(0.1,0.1,0.1,0.1))
	plot(1, xlab="", ylab="", main="", bty="n", type="n", axes=FALSE)
	if (is.null(digit))
		digit=felse(maxval>1, 0, 1)
	legValues=pretty(brks)
	legValues=legValues[which(legValues > min(brks) & legValues < max(brks))]
	shape::colorlegend(zlim=c(min(brks), max(brks)), digit=digit, col=cols, zval=legValues, posx=c(0.4, 0.5), posy = c(0.1, 0.9), xpd=TRUE)
	
	# set plotting variables
	if (!is.null(studyareaSHP)) {
		xlimit=studyareaSHP@bbox[1,]
		ylimit=studyareaSHP@bbox[2,]
	} else {
		xlimit=c(min(xmin(rast1), xmin(rast2)), max(xmax(rast1), xmax(rast2)))
		ylimit=c(min(ymin(rast1), ymin(rast2)), max(ymax(rast1), ymax(rast2)))
	}
	
	# geoplot rast 1
	par(mar=c(0.5,0.5,0.5,0.5))
	plot(1, xlim=xlimit, ylim=ylimit, type="n", axes=FALSE, bty="n")
	if (!is.null(basemapSHP)) {
		plot(basemapSHP, col="grey80", border="black", add=TRUE, lwd=10)
	}
	plotRaster(rast1, col=cols, maxpixels=maxpixels, breaks=brks, legend=FALSE, box=FALSE, axes=FALSE, ann=F, asp=1, bty="n", main=main, add=TRUE)
	if (!is.null(ptsLIST)) {Map(points, ptsLIST, col=ptsCols)}
	mtext(side=3, text=main[1], line=-1)
	
	# geoplot rast 2
	plot(1, xlim=xlimit, ylim=ylimit, type="n", axes=FALSE, bty="n")
	if (!is.null(basemapSHP)) {
		plot(basemapSHP, col="grey80", border="grey90", add=TRUE)
	}
	plotRaster(rast2, col=cols, maxpixels=maxpixels, breaks=brks, legend=FALSE, box=FALSE, axes=FALSE,  bty="n", main=main, add=TRUE)
	if (!is.null(ptsLIST)) {Map(points, ptsLIST, col=ptsCols)}
	mtext(side=3, text=main[2], line=-1)
	
	# close plotting device if saving to file
	if (!is.null(path)) {
		dev.off()
	}
}


colorlegend=function (col = femmecol(100), zlim, zlevels = 5, dz = NULL, 
    zval = NULL, zlabel = LL, log = FALSE, posx = c(0.9, 0.93), posy = c(0.05, 
        0.9), main = NULL, main.cex = 1, main.col = "black", 
    lab.col = "black", digit = 0, left = FALSE, ...) 
{
    ncol <- length(col)
    par(new = TRUE)
    omar <- nmar <- par("mar")
    nmar[c(2, 4)] <- 0
    par(mar = nmar)
    shape:::emptyplot()
    pars <- par("usr")
    dx <- pars[2] - pars[1]
    xmin <- pars[1] + posx[1] * dx
    xmax <- pars[1] + posx[2] * dx
    dy <- pars[4] - pars[3]
    ymin <- pars[3] + posy[1] * dy
    ymax <- pars[3] + posy[2] * dy
    if (!is.null(zval)) {
        zz <- zval
        dz <- NULL
    }
    if (is.null(dz) & is.null(zval)) 
        if (!is.null(zlevels)) {
            if (log) {
                zz <- 10^(pretty(log10(zlim), n = (zlevels + 
                  1)))
            }
            else zz <- pretty(zlim, n = (zlevels + 1))
        }
        else zz <- NULL
    if (!is.null(dz)) {
        if (log) 
            zz <- 10^(seq(log10(zlim[1]), log10(zlim[2]), by = dz))
        if (!log) 
            zz <- seq(zlim[1], zlim[2], by = dz)
    }
    if (log) {
        zlim <- log10(zlim)
        if (!is.null(zz)) 
            zz <- log10(zz)
    }
    zmin <- zlim[1]
    zmax <- zlim[2]
    Y <- seq(ymin, ymax, length.out = ncol + 1)
    rect(xmin, Y[-(ncol + 1)], xmax, Y[-1], col = col, border = NA)
    rect(xmin, ymin, xmax, ymax, border = lab.col)
    if (!is.null(zz)) {
        dx <- (xmax - xmin)
        dy <- (ymax - ymin)
        if (left) {
            Dx <- -dx
            pos <- 2
            xpos <- xmin + Dx * 0.5
        }
        else {
            Dx <- +dx
            pos <- 4
            xpos <- xmax + Dx * 0.5
        }
        Ypos <- ymin + (zz - zmin)/(zmax - zmin) * dy
        segments(xmin, Ypos, xmax, Ypos, col = lab.col)
        segments(xpos + Dx * 0.25, Ypos, xmin, Ypos, col = lab.col)
        if (is.null(zlabel)) {
		text(xpos, Ypos, formatC(zz, digits = digit, format = "f"), pos = pos, col = lab.col, ...)
	} else {
		text(xpos, Ypos, formatC(zlabel, digits = digit, format = "f"), pos = pos, col = lab.col, ...)
	}
    }
    if (!is.null(main)) {
        for (i in length(main):1) text(x = mean(c(xmin, xmax)), 
            y = ymax + 0.05 * (length(main) - i + 1), labels = main[i], 
            adj = c(0.5, 0.5), cex = main.cex, col = main.col)
    }
    par(new = FALSE)
    par(mar = omar)
}


### example usage
# singleMapFun(currResLYR2, studyareaSHP=extSHP, pointsSHP=currPTS, basemapSHP=basemapSHP, revpal=FALSE, digit=2, palname='Greys', letter='a', main="", res="high", colScale='quantile', path=file.path(dir(expDIR,full.names=TRUE)[i], "geoplot_cs_resistance.tiff"))