修复插值极坐标图函数以与当前 R 一起使用并（可能）使用 ggplot

Question

问题R interpolated polar contour plot 展示了一种在 R 中生成插值极坐标图的极好方法。我包括了我正在使用的稍微修改过的版本：

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
 #    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 = ceiling(max(sqrt(x^2+y^2))), 
    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 = 21, bg ="blue")
    }
    # 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)
    }
}

不幸的是，这个函数有一些错误：

a) 即使是纯粹的放射状图案，生成的图也会有失真，我不明白它的起源：

#example
r <- rep(seq(0.1, 0.9, len = 8), each = 8)
theta <- rep(seq(0, 7/4*pi, by = pi/4), times = 8)
x <- r*sin(theta)
y <- r*cos(theta)
z <- z <- rep(seq(0, 1, len = 8), each = 8)
PolarImageInterpolate(x, y, z)

为什么会在 300° 和 360° 之间摆动？ z 函数在 theta 中是不变的，所以没有理由应该有摆动。

b) 4 年后，加载的某些包已被修改，并且该功能的至少一项功能被破坏。设置interp.type = 2 应该使用薄板样条进行插值而不是基本的线性插值，但它不起作用：

> PolarImageInterpolate(x, y, z, interp.type = 2)
Warning: 
Grid searches over lambda (nugget and sill variances) with  minima at the endpoints: 
  (GCV) Generalized Cross-Validation 
   minimum at  right endpoint  lambda  =  9.493563e-06 (eff. df= 60.80002 )
predict.surface is now the function predictSurface

 Error in image.default(x = minitics, y = minitics, Mat_plot, useRaster = TRUE,  : 
  'z' must be a matrix 

第一条消息是一个警告，我并不担心，但第二条实际上是一个错误，阻止我使用薄板样条。你能帮我解决这两个问题吗？

另外，我想“升级”到使用ggplot2，所以如果你能给出一个答案，那就太好了。否则，在修复错误后，我将尝试提出一个特定问题，该问题仅要求修改函数以使其使用ggplot2。

Answer 1

对于ggplot2 解决方案，这是一个开始。 geom_raster 允许插值，但不适用于极坐标。相反，您可以使用geom_tile，不过您可能需要自己进行插值，然后再将值传递给ggplot。

重要提示：您提供的示例数据在使用geom_raster 时会出现错误，我认为这是由值的间距引起的。这是一个有效的示例集（注意，使用this blog 作为指导，尽管它现在已经过时了）：

dat_grid <-
  expand.grid(x = seq(0,350,10), y = 0:10)
dat_grid$density <- runif(nrow(dat_grid))


ggplot(dat_grid
       , aes(x = x, y = y, fill = density)) +
  geom_tile() +
  coord_polar() +
  scale_x_continuous(breaks = seq(0,360,90)) +
  scale_fill_gradient2(low = "white"
                       , mid = "yellow"
                       , high = "red3"
                       , midpoint = 0.5)

如果您使用的是原始数据，您也许可以让ggplot 为您完成这项工作。这是一个使用原始数据的示例。有很多手动修补的事情要做，但这至少是一个可选的起点：

polarData <-
  data.frame(
    theta = runif(10000, 0, 2*pi)
    , r = log(abs(rnorm(10000, 0, 10)))
  )


toCart <-
  data.frame(
    x = polarData$r * cos(polarData$theta)
    , y = polarData$r * sin(polarData$theta)
  )



axisLines <-
  data.frame(
    x = 0
    , y = 0
    , xend = max(polarData$r)*cos(seq(0, 2*pi, pi/4))
    , yend = max(polarData$r)*sin(seq(0, 2*pi, pi/4))
    , angle = paste(seq(0, 2, 1/4), "pi")  )


ticks <-
  data.frame(
    label = pretty(c(0, max(polarData$r)) )[-1]
  )


ggplot(toCart) +
  # geom_point(aes(x = x, y = y)) +
  stat_density_2d(aes(x = x, y = y
                      , fill = ..level..)
                  , geom = "polygon") +

  scale_fill_gradient(low = "white"
                      , high = "red3") +

  theme(axis.text = element_blank()
        , axis.title = element_blank()
        , axis.line = element_blank()
        , axis.ticks = element_blank()) +
  geom_segment(data = axisLines
               , aes(x = x, y = y
                     , xend = xend
                     , yend = yend)) +
  geom_label(data = axisLines
             , aes(x = xend, y = yend, label = angle)) +
  geom_label(data = ticks
             , aes(x = 0, y = label, label = label))

Answer 2

从另一篇文章中，我了解到 predict.surface 包中的函数 fields 已被弃用，它用于 PolarImageInterpolate 中的 interp.type = 2。而是引入了一个新的predictSurface函数，可以在这里使用。

例子：

r <- rep(seq(0.1, 0.9, len = 8), each = 8)
theta <- rep(seq(0, 7/4*pi, by = pi/4), times = 8)
x <- r*sin(theta)
y <- r*cos(theta)
z <- z <- rep(seq(0, 1, len = 8), each = 8)
PolarImageInterpolate(x, y, z, interp.type = 2)