Fix interpolated polar contour plot function to wo

The question R interpolated polar contour plot shows an excellent way to produce interpolated polar plots in R. I include the very slightly modified version I'm using:

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

Unfortunately, this function has a few bugs:

a) Even with a purely radial pattern, the produced plot has a distortion whose origin I don't understand:

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

why the wiggles between 300° and 360°? The z function is constant in theta, so there's no reason why there should be wiggles.

b) After 4 years, some of the packages loaded have been modified and at least one functionality of the function is broken. Setting interp.type = 2 should use thin plate splines for interpolation instead than a basic linear interpolation, but it doesn't work:

> 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

the first message is a warning and doesn't worry me, but the second one is actually an error and prevents me from using thin plate splines. Can you help me solve these two problems?

Also, I'd like to "upgrade" to using ggplot2, so if you can give an answer which does that, it would be great. Otherwise, after the bugs are fixed, I'll try asking a specific question which only asks to modify the function so that it uses ggplot2.

回答1:

For the ggplot2 solution, here is a start. geom_raster allows interpolation, but does not work for polar coordinates. Instead, you can use geom_tile, though then you may need to do the interpolation yourself before passing the values to ggplot.

Of important note: the example data you gave gives an error when working with geom_raster that I believe is caused by the spacing of the values. Here is an example set that works (note, used this blog as a guide, though it is now outdated):

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)

If you are working from raw data, you might be able to get ggplot to do the work for you. Here is an example working from raw data. There are a lot of manual tinkering things to do, but it is at least an optional starting point:

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

回答2:

From an another post, I came to know that the fucnction predict.surface from package fields is deprecated whic is used for interp.type = 2 in PolarImageInterpolate. Instead, a new predictSurface function is introduced, which can be used here.

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, interp.type = 2)