jeffreyhanson · April 21, 2015 02:46
diff --git a/geoplots.R b/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"))
	# 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"))