将等高线与等高线填充图对齐不规则网格极坐标图（半圆）

我看到有几个人回答了使用不规则网格绘图的问题。我无法使轮廓线与填充轮廓对齐。此外，需要在绘图上显示数据点位置，以及以 30 度增量显示的径向辐条，以及在 10、20 30 处显示的半圆。

Ref: 在不规则网格上绘制等高线 https://stackoverflow.com/questions/19339296/plotting-contours-on-an-irregular-grid

heading=seq(0,180,30)
speed=c(5,10,15,20,30)
mheading=matrix(heading,ncol=length(heading),nrow=length(speed),byrow=TRUE)
mspeed=matrix(speed,ncol=length(heading),nrow=length(speed),byrow=FALSE)
mag=mheading+mspeed

x=sin(mheading*pi/180)*mspeed
y=cos(mheading*pi/180)*mspeed
z=mag

library(akima)
df<-data.frame(x=x,y=y,z=z)

# interpolation
fld <- with(df, interp(x = x, y = y, z = z,
                       xo=seq(min(x),max(x),length=100),
                       yo=seq(min(y),max(y),length=100)))

filled.contour(x = fld$x,
               y = fld$y,
               z = fld$z,
               color.palette =
                 colorRampPalette(c("white", "blue")),
               xlab = "",
               ylab = "",
               main = "Max",
               key.title = title(main = "Value", cex.main = 1),
               asp=1,xlim=c(0,40),ylim=c(-30,30))

contour(x = fld$x,
               y = fld$y,
               z = fld$z,
               color.palette =
                 colorRampPalette(c("white", "blue")),
               xlab = "",
               ylab = "",
               asp=1,xlim=c(0,40),ylim=c(-30,30), add=TRUE)

按照此link https://stackoverflow.com/questions/10856882/r-interpolated-polar-contour-plot，产生以下代码/图。这是“更好”，但仍然存在问题。为什么在最小速度半径 (5) 内部存在插值数据，以及为什么轮廓填充/线没有延伸到外半径，尤其是在 90 度附近？

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?
single_point_overlay=0
outer.radius=30
spatial_res=1000 #Resolution of fitted surface
interp.type = 1
circle.rads <- pretty(c(0,outer.radius))

heading=seq(0,180,30)
speed=c(5,10,15,20,30)
mheading=matrix(heading,ncol=length(heading),nrow=length(speed),byrow=TRUE)
mspeed=matrix(speed,ncol=length(heading),nrow=length(speed),byrow=FALSE)
mag=mheading+mspeed

x=sin(mheading*pi/180)*mspeed
y=cos(mheading*pi/180)*mspeed
z=mag
extrapolate=FALSE # Should we extrapolate outside data points?

contour_levels = 8
col_levels=contour_levels
col_breaks_source=1
contour_breaks_source = 1
col = rev(heat.colors(col_levels))

minitics <- seq(-outer.radius, outer.radius, length.out = spatial_res)
xmini <- seq(min(x),max(x),length=spatial_res)
ymini <- seq(min(y),max(y),length=spatial_res)
# interpolate the data
if (interp.type ==1 ){
  #   Interp <- akima:::interp(x = x, y = y, z = z, 
  #                            extrap = extrapolate, 
  #                            xo = xmini, 
  #                            yo = ymini, 
  #                            linear = FALSE)
  #   Mat <- Interp[[3]]
  df<-data.frame(x=x,y=y,z=z)

  # interpolation
  fld <- with(df, akima:::interp(x = x, y = y, z = z,
                                 xo=xmini,
                                 yo=ymini))


  Mat_x <- fld[[1]]
  Mat_y <- fld[[2]]
  Mat_z <- fld[[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_z = tmp$z
  # mark cells outside circle as NA
  markNA <- matrix(minitics, ncol = spatial_res, nrow = spatial_res) 
  Mat_x <- markNA
  Mat_y <- t(markNA)
} else {stop("interp.type value not valid")}


# 
Mat_z[!(sqrt(Mat_x ^ 2 + Mat_y ^ 2) <= max(speed)*1.1)] <- NA 
Mat_z[!(sqrt(Mat_x ^ 2 + Mat_y ^ 2) >= min(speed))] <- NA   # <- SHOULD REMOVE INNER DATA

### 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))
  contour_breaks = seq(min(z,na.rm=TRUE),max(z,na.rm=TRUE),length.out = contour_levels+1)
}else if ((length(contour_breaks_source == 1)) & (contour_breaks_source[1] == 2)){
  contour_breaks = seq(min(Mat_z,na.rm=TRUE),max(Mat_z,na.rm=TRUE),
                       by=(max(Mat_z,na.rm=TRUE)-min(Mat_z,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_z,na.rm=TRUE),max(Mat_z,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_z
Mat_plot[which(Mat_plot<zlim[1])]=zlim[1]
Mat_plot[which(Mat_plot>zlim[2])]=zlim[2]
image(x = Mat_x, y = Mat_y, Mat_plot , 
      useRaster = TRUE, asp = 1, axes = FALSE, xlab = "", ylab = "", zlim = zlim, col = col)

# add contours if desired
if (contours){
  CL <- contourLines(x = Mat_x, y = Mat_y, Mat_z, levels = contour_breaks)
  A <- lapply(CL, function(xy){
    lines(xy$x, xy$y, col = gray(.2), lwd = .5, asp=1)
  })
}
# 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)
    # matrix(c(sin(radians), -cos(radians), cos(radians), sin(radians)), ncol = 2)
  }    

  circle <- function(x, y, rad = 1, nvert = 500, angle=360){
    rads <- seq(0,angle*pi/180,length.out = nvert)
    #     xcoords <- cos(rads) * rad + x
    #     ycoords <- sin(rads) * rad + y
    xcoords=sin(rads)*rad + x
    ycoords=cos(rads)*rad + y
    cbind(xcoords, ycoords)
  }

  # draw circles
  if (missing(circle.rads)){
    circle.rads <- pretty(c(0,outer.radius))
  }

  endAngle = 180
  for (i in circle.rads){
    lines(circle(0, 0, i, angle = endAngle), col = "#66666650")
  }

  # put on radial spoke axes:

  axis.degs <- c(0, 30, 60, 90, 120, 150)
  # axis.rads <- c(0, pi / 6, pi / 3, pi / 2, 2 * pi / 3, 5 * pi / 6)
  axis.rads <- axis.degs * pi/180
  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)){ 
    if (axis.degs[i]==0) {
      # endpoints <- zapsmall(c(RMat(axis.rads[i]) %*% matrix(c(1, 0, -1, 0) * outer.radius,ncol = 2)))
      endpoints <- zapsmall(c(RMat(axis.rads[i]) %*% matrix(c(1, 0, 0, 0) * outer.radius,ncol = 2)))
    } else if (0 < axis.degs[i] & axis.degs[i] < 90) {
      endpoints <- zapsmall(c(RMat(axis.rads[i]) %*% matrix(c(1, 0, 0, 0) * outer.radius,ncol = 2)))
    } else {
      endpoints <- zapsmall(c(RMat(axis.rads[i]) %*% matrix(c(0, 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)
    if (0 <= r.labs[i] & r.labs[i] <= 180) text(endpoints[1], endpoints[2], lab1, xpd = TRUE)
    if (0 <= l.labs[i] & l.labs[i] <= 180) 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)
  axis(2, pos = 0 * outer.radius, at = sort(union(circle.rads,-circle.rads)), labels = NA)
  text( -0.02 * outer.radius, sort(union(circle.rads, -circle.rads)),
        abs(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)
}

轮廓和颜色不对齐，因为filled.contour 生成两个图（图例和轮廓）。绘图后，这些坐标系就会丢失。 （?filled.contour）。可以通过添加相关命令来解决plot.axes争论。可以用以下方法绘制半圆draw.arc来自plotrix封装，辐条segments。最小半径内的区域可以用白色线段覆盖来表示no data.

# min distance of contours lines from center
min_dist=5

# position of spokes (degrees)
spk = seq(0,180,30)

filled.contour(x = fld$x,
               y = fld$y,
               z = fld$z,
               color.palette = colorRampPalette(c("white", "blue")),
               xlab = "",
               ylab = "",
               main = "Max",
               key.title = title(main = "Value", cex.main = 1),
               asp=1, xlim=c(0,40), ylim=c(-30,30),   frame.plot=F,
               plot.axes = {contour(fld$x, fld$y, fld$z , add=T, levels = seq(0,max(fld$z, na.rm=T),30), drawlabels=F, col=2);
                            # semi circles
                            draw.arc(x=0,y=0,radius = (1:3)*10, deg1=90, deg2=-90, col='grey');
                            # cover zone within minimum radius with (draw many closely spaced white lines
                            segments(x0 = 0, y0 = 0, x1 = sin((0:180)*pi/180)*min_dist, y1 = cos((0:180)*pi/180)*min_dist, col='white');
                            # spokes with labels
                            segments(x0 = 0, y0 = 0, x1 = sin(spk*pi/180)*30, y1 = cos(spk*pi/180)*30, col='grey');
                            text(x = sin(spk*pi/180)*30, y=cos(spk*pi/180)*30, labels = spk, pos=4, cex=0.6, xpd=NA)
                            # data points
                            points(x,y, pch=16, cex=0.6);
                            # x axis
                            axis(1);
                            # modified y axis
                            axis(2, at = axisTicks(range(y), log=F), labels = abs(axisTicks(range(y), log=F)), pos = 0);
               }
)