Ich versuche, aus interpolierten Punktdaten ein Konturpolardiagramm in R zu erstellen. Mit anderen Worten, ich habe Daten in Polarkoordinaten mit einem Größenwert, den ich zeichnen und interpolierte Werte anzeigen möchte. Ich möchte Plots wie die folgenden (in OriginPro erstellten) in Serie produzieren:

Mein nächster Versuch in R zu diesem Punkt ist im Grunde:

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

Was erzeugt so etwas wie:

Obwohl dies offensichtlich nicht das Endprodukt sein wird, ist dies der beste Weg, um Konturpolarkurven in R zu erstellen?

Ich kann zu dem Thema nur eine Archiv-Mailingliste findenDiskussion von 2008. Ich glaube, ich bin nicht ganz der Verwendung von R für die Diagramme verpflichtet (obwohl ich hier die Daten habe), aber ich bin gegen die manuelle Erstellung. Wenn es also eine andere Sprache mit dieser Fähigkeit gibt, schlagen Sie sie bitte vor (ich habe die gesehen)Python-Beispiel).

BEARBEITEN

In Bezug auf den Vorschlag mit ggplot2 - Ich kann die Routine geom_tile scheinbar nicht dazu bringen, interpolierte Daten in polar_coordinates zu zeichnen. Ich habe unten Code eingefügt, der zeigt, wo ich bin. Ich kann das Original auf kartesisch und polar zeichnen, aber ich kann nur die interpolierten Daten auf kartesisch zeichnen lassen. Ich kann die Interpolationspunkte in polar mit geom_point zeichnen, aber ich kann diesen Ansatz nicht auf geom_tile erweitern. Meine einzige Vermutung war, dass dies mit der Datenreihenfolge zusammenhängt - d. H. Geom_tile erwartet sortierte / geordnete Daten - Ich habe jede mir in den Sinn kommende Iteration ausprobiert, um die Daten ohne Änderung in aufsteigende / absteigende Azimut- und Zenitwerte zu sortieren.

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

Endgültige Ergebnisse

Ich habe schließlich den Code aus der akzeptierten Antwort (Basisgrafiken) genommen und den Code aktualisiert. Ich habe eine Spline-Interpolationsmethode für dünne Platten, eine Extrapolations- oder Nicht-Extrapolationsoption, Datenpunktüberlagerungen und die Möglichkeit hinzugefügt, kontinuierliche Farben oder segmentierte Farben für die interpolierte Oberfläche zu erstellen. Siehe die folgenden Beispiele.

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