Created
June 5, 2012 19:41
-
-
Save timriffe/2877281 to your computer and use it in GitHub Desktop.
PolarImagePlot()
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# by Tim Riffe, June 5th, 2012 | |
# arguments: | |
# Mat, a matrix of z values as follows: | |
# leftmost edge of first column = 0 degrees, rightmost edge of last column = 360 degrees | |
# columns are distributed in cells equally over the range 0 to 360 degrees, like a grid prior to transform | |
# first row is innermost circle, last row is outermost circle | |
# outer.radius, By default everything scaled to unit circle | |
# cols: color vector. default = rev(heat.colors(length(breaks)-1)) | |
# breaks: manual breaks for colors. defaults to seq(min(Mat),max(Mat),length=nbreaks) | |
# nbreaks: how many color levels are desired? | |
# axes: should circular and radial axes be drawn? radial axes are drawn at 30 degree intervals only- this could be made more flexible. | |
# circle.rads: at which radii should circles be drawn? defaults to pretty(((0:ncol(Mat)) / ncol(Mat)) * outer.radius) | |
# TODO: add color strip legend. Contours. | |
PolarImagePlot <- function(Mat, outer.radius = 1, cols, ppa = 5, breaks, nbreaks = 51, axes = TRUE, circle.rads){ | |
# the image prep | |
Mat <- Mat[,ncol(Mat):1 ] | |
radii <- ((0:ncol(Mat)) / ncol(Mat)) * outer.radius | |
# 5 points per arc will usually do | |
Npts <- ppa | |
# all the angles for which a vertex is needed | |
radians <- 2 * pi * (0:(nrow(Mat) * Npts)) / (nrow(Mat) * Npts) + pi / 2 | |
# matrix where each row is the arc corresponding to a cell | |
rad.mat <- matrix(radians[-length(radians)], ncol = Npts, byrow = TRUE)[1:nrow(Mat), ] | |
rad.mat <- cbind(rad.mat, rad.mat[c(2:nrow(rad.mat), 1), 1]) | |
# the x and y coords assuming radius of 1 | |
y0 <- sin(rad.mat) | |
x0 <- cos(rad.mat) | |
# dimension markers | |
nc <- ncol(x0) | |
nr <- nrow(x0) | |
nl <- length(radii) | |
# make a copy for each radii, redimension in sick ways | |
x1 <- aperm( x0 %o% radii, c(1, 3, 2)) | |
# the same, but coming back the other direction to close the polygon | |
x2 <- x1[, , nc:1] | |
#now stick together | |
x.array <- abind:::abind(x1[, 1:(nl - 1), ], x2[, 2:nl, ], matrix(NA, ncol = (nl - 1), nrow = nr), along = 3) | |
# final product, xcoords, is a single vector, in order, | |
# where all the x coordinates for a cell are arranged | |
# clockwise. cells are separated by NAs- allows a single call to polygon() | |
xcoords <- aperm(x.array, c(3, 1, 2)) | |
dim(xcoords) <- c(NULL) | |
# repeat for y coordinates | |
y1 <- aperm( y0 %o% radii,c(1, 3, 2)) | |
y2 <- y1[, , nc:1] | |
y.array <- abind:::abind(y1[, 1:(length(radii) - 1), ], y2[, 2:length(radii), ], matrix(NA, ncol = (length(radii) - 1), nrow = nr), along = 3) | |
ycoords <- aperm(y.array, c(3, 1, 2)) | |
dim(ycoords) <- c(NULL) | |
# sort out colors and breaks: | |
if (!missing(breaks) & !missing(cols)){ | |
if (length(breaks) - length(cols) != 1){ | |
stop("breaks must be 1 element longer than cols") | |
} | |
} | |
if (missing(breaks) & !missing(cols)){ | |
breaks <- seq(min(Mat,na.rm = TRUE), max(Mat, na.rm = TRUE), length = length(cols) + 1) | |
} | |
if (missing(cols) & !missing(breaks)){ | |
cols <- rev(heat.colors(length(breaks) - 1)) | |
} | |
if (missing(breaks) & missing(cols)){ | |
breaks <- seq(min(Mat,na.rm = TRUE), max(Mat, na.rm = TRUE), length = nbreaks) | |
cols <- rev(heat.colors(length(breaks) - 1)) | |
} | |
# get a color for each cell. Ugly, but it gets them in the right order | |
cell.cols <- as.character(cut(as.vector(Mat[nrow(Mat):1,ncol(Mat):1]), breaks = breaks, labels = cols)) | |
# start empty plot | |
plot(NULL, type = "n", ylim = c(-1, 1) * outer.radius, xlim = c(-1, 1) * outer.radius, asp = 1, axes = FALSE, xlab = "", ylab = "") | |
# draw polygons with no borders: | |
polygon(xcoords, ycoords, col = cell.cols, border = NA) | |
if (axes){ | |
# a couple internals for axis markup | |
RMat <- function(radians){ | |
matrix(c(cos(radians), sin(radians), -sin(radians), cos(radians)), ncol = 2) | |
} | |
circle <- function(x, y, rad = 1, nvert = 500){ | |
rads <- seq(0,2*pi,length.out = nvert) | |
xcoords <- cos(rads) * rad + x | |
ycoords <- sin(rads) * rad + y | |
cbind(xcoords, ycoords) | |
} | |
# draw circles | |
if (missing(circle.rads)){ | |
circle.rads <- pretty(radii) | |
} | |
for (i in circle.rads){ | |
lines(circle(0, 0, i), col = "#66666650") | |
} | |
# put on radial spoke axes: | |
axis.rads <- c(0, pi / 6, pi / 3, pi / 2, 2 * pi / 3, 5 * pi / 6) | |
r.labs <- c(90, 60, 30, 0, 330, 300) | |
l.labs <- c(270, 240, 210, 180, 150, 120) | |
for (i in 1:length(axis.rads)){ | |
endpoints <- zapsmall(c(RMat(axis.rads[i]) %*% matrix(c(1, 0, -1, 0) * outer.radius,ncol = 2))) | |
segments(endpoints[1], endpoints[2], endpoints[3], endpoints[4], col = "#66666650") | |
endpoints <- c(RMat(axis.rads[i]) %*% matrix(c(1.1, 0, -1.1, 0) * outer.radius, ncol = 2)) | |
lab1 <- bquote(.(r.labs[i]) * degree) | |
lab2 <- bquote(.(l.labs[i]) * degree) | |
text(endpoints[1], endpoints[2], lab1, xpd = TRUE) | |
text(endpoints[3], endpoints[4], lab2, xpd = TRUE) | |
} | |
axis(2, pos = -1.2 * outer.radius, at = sort(union(circle.rads,-circle.rads))) | |
} | |
invisible(list(breaks = breaks, col = cols)) | |
} | |
# example (ignore the interpolation part- do that some other way) | |
set.seed(1) | |
x <- runif(20, min = 0, max = 360) | |
y <- runif(20, min = 0, max = 40) | |
z <- rnorm(20) | |
Interp <- akima:::interp(x = x, y = y, z = z, | |
extrap = TRUE, | |
xo = seq(0, 360, length.out = 300), | |
yo = seq(0, 40, length.out = 100), | |
linear = FALSE) | |
Mat <- Interp[[3]] | |
PolarImagePlot(Mat) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Thanks Tim and Clem for the scripts and ideas. @clem: I have problems with inserting a scale.
I run the image.scale function from menugget then I wrote this into my r console `PolarImagePlot <- function(Mat,theta, outer.radius = 1, ppa = 5, cols, breaks, nbreaks = 51, axes = TRUE, circle.rads,circle.labels, title=""){
make layout to add image.scale
followed by this
layout(matrix(c(1,2), nrow=1, ncol=2),width=c(6,1), height=c(4,4))then the function from Tim
# the image prepMat <- Mat[,ncol(Mat):1 ]
radii <- ((0:ncol(Mat)) / ncol(Mat)) * outer.radius
par(mar=c(1,1,1,5))
image.scale(Mat, col=cols, breaks=breaks, horiz=FALSE, yaxt="n", zlim=range(Mat, na.rm=TRUE))
axis(4,at=round(breaks[seq(1,length(breaks),length=10)], digits=2), las=2)
box()
}
When I run this function with the provided example data from Tim
# example (ignore the interpolation part- do that some other way)set.seed(1)
x <- runif(20, min = 0, max = 360)
y <- runif(20, min = 0, max = 40)
z <- rnorm(20)
Interp <- akima:::interp(x = x, y = y, z = z,
extrap = TRUE,
xo = seq(0, 360, length.out = 300),
yo = seq(0, 40, length.out = 100),
linear = FALSE)
Mat <- Interp[[3]]
PolarImagePlot(Mat)
I got the following Errormassage:
Error in .External.graphics(C_layout, num.rows, num.cols, mat, as.integer(num.figures), :invalid graphics state`
I am a R beginner and wnat to plot bidirectional data. Can you please give me hint what I am doing wrong in applying your code? sessioniInfo()