R 内插极坐标等值线图 [英] R interpolated polar contour plot
本文介绍了R 内插极坐标等值线图的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在尝试根据内插点数据编写 R 中的等高线极坐标图.换句话说，我有极坐标中的数据，我想绘制和显示插值的幅度值.我想批量生产类似于以下的图(在 OriginPro 中生产):
我在 R 中最接近的尝试基本上是:
### 转换极坐标 ->大车# 去做 #### 虚拟数据x = 范数(20)y = 范数(20)z = 范数(20)### 插值图书馆(akima)tmp = interp(x,y,z)### 绘图插值图书馆(领域)图像.plot(tmp)### 去做 ####关闭所有轴#在顶部绘制极轴产生如下结果:
虽然这显然不是最终产品，但这是在 R 中创建等高线极坐标图的最佳方式吗?
除了存档邮件列表之外，我找不到关于该主题的任何内容
[[我原来的回答如下]]
我认为你的问题对我自己很有教育意义，所以我接受了挑战并提出了以下不完整的函数.它的工作原理与 image() 类似，需要一个矩阵作为其主要输入，然后吐出类似于上面示例的内容，减去轮廓线.[[我在 6 月 6 日编辑代码后注意到它没有按照我声称的顺序进行绘制.固定的.目前正在研究等高线和图例.]]
# 参数:# Mat，z 值矩阵如下:# 第一列的最左边 = 0 度，最后一列的最右边 = 360 度# 列在 0 到 360 度范围内均匀分布在单元格中，就像变换之前的网格一样# 第一行是最内圈，最后一行是最外圈# external.radius, 默认情况下一切都缩放到单位圆# ppa:每条弧每个单元格的点数.如果您的矩阵很小，请将其放大以获得漂亮的曲线# cols:颜色向量.默认 = rev(heat.colors(length(breaks)-1))# 中断:颜色的手动中断.默认为 seq(min(Mat),max(Mat),length=nbreaks)# nbreaks:需要多少颜色级别?# 轴:是否应该绘制圆形和径向轴?径向轴仅以 30 度的间隔绘制 - 这可以变得更加灵活.# circle.rads:圆应该在哪个半径上画?默认为漂亮(((0:ncol(Mat))/ncol(Mat))*外层半径)# TODO:添加颜色条图例.PolarImagePlot <- function(Mat，outer.radius = 1，ppa = 5，cols，breaks，nbreaks = 51，axes = TRUE，circle.rads){# 图像准备Mat <- Mat[, ncol(Mat):1]半径 <- ((0:ncol(Mat))/ncol(Mat)) * outer.radius# 每弧 5 点通常就可以了Npts <- ppa# 需要顶点的所有角度弧度 <- 2 * pi * (0:(nrow(Mat) * Npts))/(nrow(Mat) * Npts) + pi/2# 矩阵，其中每一行是对应于一个单元格的弧rad.mat <-矩阵(弧度[-长度(弧度)]，ncol = Npts，byrow = TRUE)[1:nrow(Mat)，]rad.mat <- cbind(rad.mat, rad.mat[c(2:nrow(rad.mat), 1), 1])# x 和 y 坐标假设半径为 1y0 <- sin(rad.mat)x0 <- cos(rad.mat)# 维度标记nc <- ncol(x0)nr <- nrow(x0)nl <- 长度(半径)# 为每个半径制作一个副本，以病态的方式重新标注x1 <- aperm( x0 %o% 半径，c(1, 3, 2))# 相同，但从另一个方向返回以关闭多边形x2 <- x1[, , nc:1]#现在团结起来x.array <- abind:::abind(x1[, 1:(nl - 1), ], x2[, 2:nl, ], matrix(NA, ncol = (nl - 1), nrow = nr), 沿 = 3)# 最终产品 xcoords 是单个向量，按顺序，# 一个单元格的所有 x 坐标的排列位置# 顺时针.单元格由 NAs 分隔 - 允许一次调用多边形()xcoords <- aperm(x.array, c(3, 1, 2))暗淡(xcoords) <- c(NULL)# 重复 y 坐标y1 <- aperm( y0 %o% 半径,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), 沿 = 3)ycoords <- aperm(y.array, c(3, 1, 2))暗淡(ycoords) <- c(NULL)# 整理颜色和休息:如果 (!missing(breaks) & !missing(cols)){if (length(breaks) - length(cols) != 1){stop("breaks 必须比 cols 长 1 个元素")}}如果(缺失(中断)& ！missing(cols)){中断 <- seq(min(Mat,na.rm = TRUE), max(Mat, na.rm = TRUE), length = length(cols) + 1)}如果(缺少(列)& ！缺少(中断)){cols <- rev(heat.colors(length(breaks) - 1))}如果(缺失(中断)和缺失(列)){中断 <- seq(min(Mat,na.rm = TRUE), max(Mat, na.rm = TRUE), length = nbreaks)cols <- rev(heat.colors(length(breaks) - 1))}# 获取每个单元格的颜色.丑陋，但它使它们以正确的顺序排列cell.cols <- as.character(cut(as.vector(Mat[nrow(Mat):1,ncol(Mat):1])，break = break，labels = cols))# 开始空图plot(NULL, type = "n", ylim = c(-1, 1) * outer.radius, xlim = c(-1, 1) * outer.radius, asp = 1, axes = FALSE, xlab = "", ylab = "")# 绘制没有边界的多边形:多边形(xcoords，ycoords，col = cell.cols，border = NA)如果(轴){# 轴标记的几个内部结构.RMat <- 函数(弧度){矩阵(c(cos(弧度)，sin(弧度)，-sin(弧度)，cos(弧度))，ncol = 2)}圆<-函数(x，y，rad = 1，nvert = 500){rads <- seq(0,2*pi,length.out = nvert)xcoords <- cos(rads) * rad + xycoords <- sin(rads) * rad + ycbind(xcoords，ycoords)}# 画圆圈如果(缺少(circle.rads)){circle.rads <- 漂亮(半径)}for (i in circle.rads){线(圆圈(0、0、i)、col =#66666650")}# 放在径向辐条轴上: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)){端点 <- zapsmall(c(RMat(axis.rads[i]) %*% matrix(c(1, 0, -1, 0) * outer.radius,ncol = 2)))段(端点[1]，端点[2]，端点[3]，端点[4]，col =#66666650")端点 <- 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)文本(端点[1]，端点[2]，lab1，xpd = TRUE)文本(端点[3]，端点[4]，lab2，xpd = TRUE)}轴(2，pos = -1.2 *外层半径，at = sort(union(circle.rads，-circle.rads)))}隐形(列表(中断 = 中断，col = cols))}我不知道如何在极地表面上正确插值，因此假设您可以实现这一点并将数据放入矩阵中，那么此函数将为您绘制它.每个单元格都被绘制，就像 image() 一样，但内部的单元格很小.举个例子:
set.seed(1)x <- runif(20, min = 0, max = 360)y <- runif(20, min = 0, max = 40)z <- 范数(20)Interp <- akima:::interp(x = x, y = y, z = z,额外=真，xo = seq(0, 360, length.out = 300),yo = seq(0, 40, length.out = 100),线性 = 假)Mat <- Interp[[3]]PolarImagePlot(Mat)无论如何，请随意修改它并随心所欲地使用它.代码在 Github 上可用:https://gist.github.com/2877281
I'm attempting to script a contour polar plot in R from interpolated point data. In other words, I have data in polar coordinates with a magnitude value I would like to plot and show interpolated values. I'd like to mass produce plots similar to the following (produced in OriginPro):
My closest attempt in R to this point is basically:
### Convert polar -> cart
# ToDo #
### Dummy data
x = rnorm(20)
y = rnorm(20)
z = rnorm(20)
### Interpolate
library(akima)
tmp = interp(x,y,z)
### Plot interpolation
library(fields)
image.plot(tmp)
### ToDo ###
#Turn off all axis
#Plot polar axis ontop
Which produces something like:
While this is obviously not going to be the final product, is this the best way to go about creating contour polar plots in R?
I can't find anything on the topic other than an archive mailing list dicussion from 2008. I guess I'm not fully dedicated to using R for the plots (though that is where I have the data), but I am opposed to manual creation. So, if there is another language with this capability, please suggest it (I did see the Python example).
EDIT
Regarding the suggestion using ggplot2 - I can't seem to get the geom_tile routine to plot interpolated data in polar_coordinates. I have included code below that illustrates where I'm at. I can plot the original in Cartesian and polar, but I can only get the interpolated data to plot in Cartesian. I can plot the interpolation points in polar using geom_point, but I can't extend that approach to geom_tile. My only guess was that this is related to data order - i.e. geom_tile is expecting sorted/ordered data - I've tried every iteration I can think of sorting the data into ascending/descending azimuth and zenith with no change.
## Libs
library(akima)
library(ggplot2)
## Sample data in az/el(zenith)
tmp = seq(5,355,by=10)
geoms <- data.frame(az = tmp,
zen = runif(length(tmp)),
value = runif(length(tmp)))
geoms$az_rad = geoms$az*pi/180
## These points plot fine
ggplot(geoms)+geom_point(aes(az,zen,colour=value))+
coord_polar()+
scale_x_continuous(breaks=c(0,45,90,135,180,225,270,315,360),limits=c(0,360))+
scale_colour_gradient(breaks=seq(0,1,by=.1),low="black",high="white")
## Need to interpolate - most easily done in cartesian
x = geoms$zen*sin(geoms$az_rad)
y = geoms$zen*cos(geoms$az_rad)
df.ptsc = data.frame(x=x,y=y,z=geoms$value)
intc = interp(x,y,geoms$value,
xo=seq(min(x), max(x), length = 100),
yo=seq(min(y), max(y), length = 100),linear=FALSE)
df.intc = data.frame(expand.grid(x=intc$x,y=intc$y),
z=c(intc$z),value=cut((intc$z),breaks=seq(0,1,.1)))
## This plots fine in cartesian coords
ggplot(df.intc)+scale_x_continuous(limits=c(-1.1,1.1))+
scale_y_continuous(limits=c(-1.1,1.1))+
geom_point(data=df.ptsc,aes(x,y,colour=z))+
scale_colour_gradient(breaks=seq(0,1,by=.1),low="white",high="red")
ggplot(df.intc)+geom_tile(aes(x,y,fill=z))+
scale_x_continuous(limits=c(-1.1,1.1))+
scale_y_continuous(limits=c(-1.1,1.1))+
geom_point(data=df.ptsc,aes(x,y,colour=z))+
scale_colour_gradient(breaks=seq(0,1,by=.1),low="white",high="red")
## Convert back to polar
int_az = atan2(df.intc$x,df.intc$y)
int_az = int_az*180/pi
int_az = unlist(lapply(int_az,function(x){if(x<0){x+360}else{x}}))
int_zen = sqrt(df.intc$x^2+df.intc$y^2)
df.intp = data.frame(az=int_az,zen=int_zen,z=df.intc$z,value=df.intc$value)
## Just to check
az = atan2(x,y)
az = az*180/pi
az = unlist(lapply(az,function(x){if(x<0){x+360}else{x}}))
zen = sqrt(x^2+y^2)
## The conversion looks correct [[az = geoms$az, zen = geoms$zen]]
## This plots the interpolated locations
ggplot(df.intp)+geom_point(aes(az,zen))+coord_polar()
## This doesn't track to geom_tile
ggplot(df.intp)+geom_tile(aes(az,zen,fill=value))+coord_polar()
Final Results
I finally took code from the accepted answer (base graphics) and updated the code. I added a thin plate spline interpolation method, an option to extrapolate or not, data point overlays, and the ability to do continuous colors or segmented colors for the interpolated surface. See the examples below.
PolarImageInterpolate <- function(
### Plotting data (in cartesian) - will be converted to polar space.
x, y, z,
### Plot component flags
contours=TRUE, # Add contours to the plotted surface
legend=TRUE, # Plot a surface data legend?
axes=TRUE, # Plot axes?
points=TRUE, # Plot individual data points
extrapolate=FALSE, # Should we extrapolate outside data points?
### Data splitting params for color scale and contours
col_breaks_source = 1, # Where to calculate the color brakes from (1=data,2=surface)
# If you know the levels, input directly (i.e. c(0,1))
col_levels = 10, # Number of color levels to use - must match length(col) if
#col specified separately
col = rev(heat.colors(col_levels)), # Colors to plot
contour_breaks_source = 1, # 1=z data, 2=calculated surface data
# If you know the levels, input directly (i.e. c(0,1))
contour_levels = col_levels+1, # One more contour break than col_levels (must be
# specified correctly if done manually
### Plotting params
outer.radius = round_any(max(sqrt(x^2+y^2)),5,f=ceiling),
circle.rads = pretty(c(0,outer.radius)), #Radius lines
spatial_res=1000, #Resolution of fitted surface
single_point_overlay=0, #Overlay "key" data point with square
#(0 = No, Other = number of pt)
### Fitting parameters
interp.type = 1, #1 = linear, 2 = Thin plate spline
lambda=0){ #Used only when interp.type = 2
minitics <- seq(-outer.radius, outer.radius, length.out = spatial_res)
# interpolate the data
if (interp.type ==1 ){
Interp <- akima:::interp(x = x, y = y, z = z,
extrap = extrapolate,
xo = minitics,
yo = minitics,
linear = FALSE)
Mat <- Interp[[3]]
}
else if (interp.type == 2){
library(fields)
grid.list = list(x=minitics,y=minitics)
t = Tps(cbind(x,y),z,lambda=lambda)
tmp = predict.surface(t,grid.list,extrap=extrapolate)
Mat = tmp$z
}
else {stop("interp.type value not valid")}
# mark cells outside circle as NA
markNA <- matrix(minitics, ncol = spatial_res, nrow = spatial_res)
Mat[!sqrt(markNA ^ 2 + t(markNA) ^ 2) < outer.radius] <- NA
### Set contour_breaks based on requested source
if ((length(contour_breaks_source == 1)) & (contour_breaks_source[1] == 1)){
contour_breaks = seq(min(z,na.rm=TRUE),max(z,na.rm=TRUE),
by=(max(z,na.rm=TRUE)-min(z,na.rm=TRUE))/(contour_levels-1))
}
else if ((length(contour_breaks_source == 1)) & (contour_breaks_source[1] == 2)){
contour_breaks = seq(min(Mat,na.rm=TRUE),max(Mat,na.rm=TRUE),
by=(max(Mat,na.rm=TRUE)-min(Mat,na.rm=TRUE))/(contour_levels-1))
}
else if ((length(contour_breaks_source) == 2) & (is.numeric(contour_breaks_source))){
contour_breaks = pretty(contour_breaks_source,n=contour_levels)
contour_breaks = seq(contour_breaks_source[1],contour_breaks_source[2],
by=(contour_breaks_source[2]-contour_breaks_source[1])/(contour_levels-1))
}
else {stop("Invalid selection for "contour_breaks_source"")}
### Set color breaks based on requested source
if ((length(col_breaks_source) == 1) & (col_breaks_source[1] == 1))
{zlim=c(min(z,na.rm=TRUE),max(z,na.rm=TRUE))}
else if ((length(col_breaks_source) == 1) & (col_breaks_source[1] == 2))
{zlim=c(min(Mat,na.rm=TRUE),max(Mat,na.rm=TRUE))}
else if ((length(col_breaks_source) == 2) & (is.numeric(col_breaks_source)))
{zlim=col_breaks_source}
else {stop("Invalid selection for "col_breaks_source"")}
# begin plot
Mat_plot = Mat
Mat_plot[which(Mat_plot<zlim[1])]=zlim[1]
Mat_plot[which(Mat_plot>zlim[2])]=zlim[2]
image(x = minitics, y = minitics, Mat_plot , useRaster = TRUE, asp = 1, axes = FALSE, xlab = "", ylab = "", zlim = zlim, col = col)
# add contours if desired
if (contours){
CL <- contourLines(x = minitics, y = minitics, Mat, levels = contour_breaks)
A <- lapply(CL, function(xy){
lines(xy$x, xy$y, col = gray(.2), lwd = .5)
})
}
# add interpolated point if desired
if (points){
points(x,y,pch=4)
}
# add overlay point (used for trained image marking) if desired
if (single_point_overlay!=0){
points(x[single_point_overlay],y[single_point_overlay],pch=0)
}
# add radial axes if desired
if (axes){
# 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(c(0,outer.radius))
}
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.25 * outer.radius, at = sort(union(circle.rads,-circle.rads)), labels = NA)
text( -1.26 * outer.radius, sort(union(circle.rads, -circle.rads)),sort(union(circle.rads, -circle.rads)), xpd = TRUE, pos = 2)
}
# add legend if desired
# this could be sloppy if there are lots of breaks, and that's why it's optional.
# another option would be to use fields:::image.plot(), using only the legend.
# There's an example for how to do so in its documentation
if (legend){
library(fields)
image.plot(legend.only=TRUE, smallplot=c(.78,.82,.1,.8), col=col, zlim=zlim)
# ylevs <- seq(-outer.radius, outer.radius, length = contour_levels+ 1)
# #ylevs <- seq(-outer.radius, outer.radius, length = length(contour_breaks))
# rect(1.2 * outer.radius, ylevs[1:(length(ylevs) - 1)], 1.3 * outer.radius, ylevs[2:length(ylevs)], col = col, border = NA, xpd = TRUE)
# rect(1.2 * outer.radius, min(ylevs), 1.3 * outer.radius, max(ylevs), border = "#66666650", xpd = TRUE)
# text(1.3 * outer.radius, ylevs[seq(1,length(ylevs),length.out=length(contour_breaks))],round(contour_breaks, 1), pos = 4, xpd = TRUE)
}
}
解决方案
[[major edit]] I was finally able to add contour lines to my original attempt, but since the two sides of the original matrix that gets contorted don't actually touch, the lines don't match up between 360 and 0 degree. So I've totally rethought the problem, but leave the original post below because it was still kind of cool to plot a matrix that way. The function I'm posting now takes x,y,z and several optional arguments, and spits back something pretty darn similar to your desired examples, radial axes, legend, contour lines and all:
PolarImageInterpolate <- function(x, y, z, outer.radius = 1,
breaks, col, nlevels = 20, contours = TRUE, legend = TRUE,
axes = TRUE, circle.rads = pretty(c(0,outer.radius))){
minitics <- seq(-outer.radius, outer.radius, length.out = 1000)
# interpolate the data
Interp <- akima:::interp(x = x, y = y, z = z,
extrap = TRUE,
xo = minitics,
yo = minitics,
linear = FALSE)
Mat <- Interp[[3]]
# mark cells outside circle as NA
markNA <- matrix(minitics, ncol = 1000, nrow = 1000)
Mat[!sqrt(markNA ^ 2 + t(markNA) ^ 2) < outer.radius] <- NA
# sort out colors and breaks:
if (!missing(breaks) & !missing(col)){
if (length(breaks) - length(col) != 1){
stop("breaks must be 1 element longer than cols")
}
}
if (missing(breaks) & !missing(col)){
breaks <- seq(min(Mat,na.rm = TRUE), max(Mat, na.rm = TRUE), length = length(col) + 1)
nlevels <- length(breaks) - 1
}
if (missing(col) & !missing(breaks)){
col <- rev(heat.colors(length(breaks) - 1))
nlevels <- length(breaks) - 1
}
if (missing(breaks) & missing(col)){
breaks <- seq(min(Mat,na.rm = TRUE), max(Mat, na.rm = TRUE), length = nlevels + 1)
col <- rev(heat.colors(nlevels))
}
# if legend desired, it goes on the right and some space is needed
if (legend) {
par(mai = c(1,1,1.5,1.5))
}
# begin plot
image(x = minitics, y = minitics, t(Mat), useRaster = TRUE, asp = 1,
axes = FALSE, xlab = "", ylab = "", col = col, breaks = breaks)
# add contours if desired
if (contours){
CL <- contourLines(x = minitics, y = minitics, t(Mat), levels = breaks)
A <- lapply(CL, function(xy){
lines(xy$x, xy$y, col = gray(.2), lwd = .5)
})
}
# add radial axes if desired
if (axes){
# 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(c(0,outer.radius))
}
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)), labels = NA)
text( -1.21 * outer.radius, sort(union(circle.rads, -circle.rads)),sort(union(circle.rads, -circle.rads)), xpd = TRUE, pos = 2)
}
# add legend if desired
# this could be sloppy if there are lots of breaks, and that's why it's optional.
# another option would be to use fields:::image.plot(), using only the legend.
# There's an example for how to do so in its documentation
if (legend){
ylevs <- seq(-outer.radius, outer.radius, length = nlevels + 1)
rect(1.2 * outer.radius, ylevs[1:(length(ylevs) - 1)], 1.3 * outer.radius, ylevs[2:length(ylevs)], col = col, border = NA, xpd = TRUE)
rect(1.2 * outer.radius, min(ylevs), 1.3 * outer.radius, max(ylevs), border = "#66666650", xpd = TRUE)
text(1.3 * outer.radius, ylevs,round(breaks, 1), pos = 4, xpd = TRUE)
}
}
# Example
set.seed(10)
x <- rnorm(20)
y <- rnorm(20)
z <- rnorm(20)
PolarImageInterpolate(x,y,z, breaks = seq(-2,8,by = 1))
code available here: https://gist.github.com/2893780
[[my original answer follows]]
I thought your question would be educational for myself, so I took up the challenge and came up with the following incomplete function. It works similar to image(), wants a matrix as its primary input, and spits back something similar to your example above, minus the contour lines.
[[I edited the code June 6th after noticing that it didn't plot in the order I claimed. Fixed. Currently working on contour lines and legend.]]
# 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
# ppa: points per cell per arc. If your matrix is little, make it larger for a nice curve
# 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.
PolarImagePlot <- function(Mat, outer.radius = 1, ppa = 5, cols, 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))
}
I don't know how to interpolate properly over a polar surface, so assuming you can achieve that and get your data into a matrix, then this function will get it plotted for you. Each cell is drawn, as with image(), but the interior ones are teeny tiny. Here's an example:
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)
By all means, feel free to modify this and do with it what you will. Code is available on Github here: https://gist.github.com/2877281
这篇关于R 内插极坐标等值线图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
