studio - superponer graficas en r ggplot

R diagrama de contorno polar interpolado

Question

Intento crear un guión de un diagrama polar de contorno en R a partir de datos de puntos interpolados. En otras palabras, tengo datos en coordenadas polares con un valor de magnitud que me gustaría trazar y mostrar valores interpolados. Me gustaría producir parcelas similares a las siguientes (producidas en OriginPro):

Mi intento más cercano en R hasta este punto es básicamente:

### 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 produce algo como:

Si bien esto obviamente no va a ser el producto final, ¿es esta la mejor manera de crear trazados contornos polares en R?

No puedo encontrar nada sobre el tema que no sea una discusión de la lista de correo de archivo desde 2008 . Supongo que no estoy totalmente dedicado a usar R para las tramas (aunque ahí es donde tengo los datos), pero me opongo a la creación manual. Por lo tanto, si hay otro idioma con esta capacidad, sugiéralo (vi el ejemplo de Python ).

EDITAR

En cuanto a la sugerencia de usar ggplot2, parece que no puedo obtener la rutina geom_tile para trazar los datos interpolados en polar_coordinates. He incluido el código a continuación que ilustra dónde estoy. Puedo trazar el original en cartesiano y polar, pero solo puedo obtener los datos interpolados para trazarlos en cartesiano. Puedo trazar los puntos de interpolación en polar usando geom_point, pero no puedo extender ese enfoque a geom_tile. Mi única conjetura es que esto está relacionado con el orden de los datos, es decir, geom_tile está esperando datos ordenados / ordenados. He intentado cada iteración en la que puedo pensar clasificando los datos en azimut ascendente / descendente y cenit sin cambios.

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

Finalmente tomé el código de la respuesta aceptada (gráficos básicos) y actualicé el código. Agregué un método de interpolación de spline de placa delgada, una opción para extrapolar o no, superposiciones de puntos de datos, y la capacidad de hacer colores continuos o colores segmentados para la superficie interpolada. Vea los ejemplos a continuación.

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

Answer 1

Diagrama objetivo

Código de ejemplo

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 tiene muchas opciones para explorar escalas de color, anotaciones, etc., pero esto debería comenzar.

El crédito a esta respuesta de Andrie de Vries me llevó a esta solución.

Answer 2

[[edición principal]] Finalmente pude agregar líneas de contorno a mi intento original, pero como los dos lados de la matriz original que se contorsionan no se tocan, las líneas no coinciden entre 360 ​​y 0 grados. Así que he reconsiderado por completo el problema, pero dejo la publicación original a continuación porque todavía era genial trazar una matriz de esa manera. La función que estoy publicando ahora toma x, y, z y varios argumentos opcionales, y escupe algo bastante similar a los ejemplos que deseas, ejes radiales, leyendas, líneas de contorno y todo:

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

código disponible aquí: https://gist.github.com/2893780

[[mi respuesta original sigue]]

Pensé que tu pregunta sería educativa para mí, así que acepté el desafío y se me ocurrió la siguiente función incompleta. Funciona de manera similar a image() , quiere una matriz como su entrada principal, y devuelve algo similar a su ejemplo anterior, menos las líneas de contorno. [[Edité el código el 6 de junio después de darme cuenta de que no tramaba en el orden que afirmaba. Fijo. Actualmente trabajando en curvas de nivel y leyenda.]]

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

No sé cómo interpolar correctamente sobre una superficie polar, por lo que suponiendo que pueda lograr eso y obtener sus datos en una matriz, entonces esta función lo trazará para usted. Cada celda se dibuja, como con la image() , pero las interiores son minúsculas. Aquí hay un ejemplo:

    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)

Por supuesto, siéntete libre de modificar esto y hacer con él lo que quieras. El código está disponible en Github aquí: https://gist.github.com/2877281