Eu estou tentando rotear um gráfico polar de contorno em R a partir de dados de pontos interpolados. Em outras palavras, tenho dados em coordenadas polares com um valor de magnitude que gostaria de plotar e mostrar valores interpolados. Eu gostaria de produzir em massa lotes semelhantes aos seguintes (produzidos no OriginPro):

Minha tentativa mais próxima em R para este ponto é basicamente:

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

Que produz algo como:

Embora, obviamente, este não seja o produto final, essa é a melhor maneira de criar contornos polares de contorno em R?

Eu não consigo encontrar nada sobre o tema que não seja uma lista de discussão de arquivodiscussão de 2008. Eu acho que não estou totalmente dedicado ao uso de R para os gráficos (embora seja onde eu tenho os dados), mas sou contra a criação manual. Então, se houver outro idioma com essa capacidade, por favor sugira (eu vi oExemplo do Python).

EDITAR

Em relação à sugestão usando ggplot2 - não consigo fazer com que a rotina geom_tile plote dados interpolados em polar_coordinates. Eu incluí o código abaixo que ilustra onde eu estou. Eu posso traçar o original em cartesiano e polar, mas só posso obter os dados interpolados para traçar em cartesiano. Eu posso plotar os pontos de interpolação em polar usando o ponto_gal, mas não posso estender essa abordagem para o geom_tile. Meu único palpite é que isso está relacionado à ordem de dados - ou seja, geom_tile está esperando dados ordenados / ordenados - tentei todas as iterações que posso pensar em classificar os dados em ascendente / decrescente azimute e zenith sem alteração.

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

Resultados finais

Eu finalmente peguei o código da resposta aceita (gráficos de base) e atualizei o código. Eu adicionei um método de interpolação de spline de placa fina, uma opção para extrapolar ou não, sobreposições de pontos de dados e a capacidade de fazer cores contínuas ou segmentadas para a superfície interpolada. Veja os exemplos abaixo.

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