Skip to content

Instantly share code, notes, and snippets.

@jeffreyhanson

jeffreyhanson/geoplots.R

Created April 21, 2015 02:46
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save jeffreyhanson/b2c16ad6922fb7c7ba0d to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save jeffreyhanson/b2c16ad6922fb7c7ba0d to your computer and use it in GitHub Desktop.
Download ZIP
generic geoplotting functions
Raw
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"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment