-
-
Save timriffe/2877281 to your computer and use it in GitHub Desktop.
# 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) |
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.scalefollowed 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 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))` and finally the last part of your comment `#Add scale
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()
R version 3.4.0 Patched (2017-06-13 r72789)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
Thanks a lot for sharing your script! This was very helpful for our bidirectional reflectance project!!!
I used your R script and added an Color Scale using the image.scale function from http://menugget.blogspot.de/2011/08/adding-scale-to-image-plot.html within the PolarImagePLot function:
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
layout(matrix(c(1,2), nrow=1, ncol=2),width=c(6,1), height=c(4,4))
first image
par(mar=c(1,1,2,1))
... your function
#Add scale
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()
}
Regards,
Clem