R interpolées polaire tracé de contour

Que j'essaie de faire un script de contour, tracé polaire dans la R de interpolées point de données. En d'autres termes, j'ai des données en coordonnées polaires avec une grandeur de valeur, je voudrais de l'intrigue et de montrer les valeurs interpolées. J'aimerais produire en masse des parcelles semblable à la suivante (réalisée en OriginPro):

R interpolées polaire tracé de contour

Mes plus proches tentent de R pour ce point, c'est en gros:

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

Qui produit quelque chose comme:
R interpolées polaire tracé de contour

Alors que ce n'est évidemment pas le produit final, est-ce la meilleure façon d'aller sur la création de contour des tracés polaires dans R?

Je ne trouve rien sur le sujet autre qu'une archive de la liste de diffusion discussions à partir de 2008. Je suppose que je ne suis pas entièrement dédié à l'utilisation de R pour les parcelles (même si c'est là que j'ai les données), mais je suis opposé à la création manuelle. Donc, si il ya une autre langue grâce à cette fonctionnalité, s'il vous plaît suggérer qu'il (je l'ai fait voir la Exemple Python).

MODIFIER

S'agissant de la suggestion à l'aide de ggplot2 - je n'arrive pas à obtenir le geom_tile de routine pour le tracé de données interpolées dans polar_coordinates. J'ai inclus le code ci-dessous qui illustre où je suis. Je peux tracer l'origine en coordonnées Cartésiennes et polaires, mais je ne peux faire l'interpolation des données de tracé en coordonnées Cartésiennes. Je peux tracer l'interpolation des points en polaire à l'aide de geom_point, mais je ne peux pas étendre cette approche à geom_tile. Mon seul suppose que cela est lié à des données de commande (geom_tile s'attend à trier/a ordonné de données, j'ai essayé à chaque itération, je pense, de trier les données en ordre croissant/décroissant d'azimut et d'zénith aucun changement.

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

Résultats Définitifs

Finalement j'ai pris le code de la accepté de répondre (de graphismes de base) et de la mise à jour du code. J'ai ajouté une fine plaque d'interpolation spline méthode, une option permettant d'extrapoler ou pas, le point de données des superpositions, et la capacité de le faire en continu couleurs ou segmentées couleurs de la surface interpolée. Voir les exemples ci-dessous.

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

R interpolées polaire tracé de contour
R interpolées polaire tracé de contour
R interpolées polaire tracé de contour

  • génial, merci pour la publication de votre résultat final
InformationsquelleAutor Nate | 2012-06-01