R插值极坐标轮廓图

2 回答

• 5

  [[主编辑]]我终于能够将轮廓线添加到我原来的尝试中，但由于原始矩阵的两边被扭曲实际上没有触及，所以这些线在360度和0度之间不匹配 . 所以我完全重新思考了这个问题，但是保留下面的原始帖子，因为以这种方式绘制矩阵仍然很酷 . 我现在发布的函数采用x，y，z和几个可选参数，然后吐出一些非常类似于你想要的例子，径向轴，图例，轮廓线和所有：

  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))
  

  代码在这里：https://gist.github.com/2893780

  

  [[我原来的答案如下]]

  我认为你的问题对我自己很有教育意义，所以我接受了挑战并提出了以下不完整的功能 . 它与 image() 的工作方式类似，想要一个矩阵作为它的主要输入，然后回吐类似于上面例子的东西，减去轮廓线 . [[我在注意到它没有按照我声称的顺序绘制后，于6月6日编辑了代码 . 固定 . 目前正在研究轮廓线和图例 . ]]

  # 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))
      }
  

  我不知道如何在极性表面上正确插值，因此假设您可以实现并将数据转换为矩阵，那么此函数将为您绘制 . 绘制每个单元格，就像 image() 一样，但是内部单元格很小 . 这是一个例子：

  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)
  

  

  无论如何，随意修改它并随意使用它 . 代码在Github上提供：https://gist.github.com/2877281

  回复于 2024-04-29T02:22:06+08:00
• 11

  目标图

  

  示例代码

  library(akima) 
  library(ggplot2) 
  
  x = rnorm(20)
  y = rnorm(20)
  z = rnorm(20)
  
  t. = interp(x,y,z)
  t.df <- data.frame(t.)
  
  gt <- data.frame( expand.grid(X1=t.$x, 
                                X2=t.$y), 
                    z=c(t.$z), 
                    value=cut(c(t.$z), 
                              breaks=seq(-1,1,0.25)))
  
  p <- ggplot(gt) + 
      geom_tile(aes(X1,X2,fill=value)) + 
      geom_contour(aes(x=X1,y=X2,z=z), colour="black") + 
      coord_polar()
  p <- p + scale_fill_brewer()
  p
  

  ggplot2 然后有很多选项来探索重新色标，注释等，但这应该让你开始 .

  感谢this answer by Andrie de Vries，这让我得到了这个解决方案 .

  回复于 2024-04-29T02:22:06+08:00