title: R interpolated polar contour plot
body: <p>I'm attempting to script a contour polar plot in R from interpolated point data.  In other words, I have data in polar coordinates with a magnitude value I would like to plot and show interpolated values.  I'd like to mass produce plots similar to the following (produced in OriginPro):</p>
  
  <p><img src="https://i.stack.imgur.com/L7l7x.png" alt="OriginPro contour polar plot"></p>
  
  <p>My closest attempt in R to this point is basically:</p>
  
  <pre><code>### Convert polar -&gt; 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
  </code></pre>
  
  <p>Which produces something like:
  <img src="https://i.stack.imgur.com/m9gG6.png" alt="Dummy R interpolated plot"></p>
  
  <p>While this is obviously not going to be the final product, is this the best way to go about creating contour polar plots in R?  </p>
  
  <p>I can't find anything on the topic other than an archive mailing list <a href="https://stat.ethz.ch/pipermail/r-help/2008-December/182166.html" rel="noreferrer">dicussion from 2008</a>.  I guess I'm not fully dedicated to using R for the plots (though that is where I have the data), but I am opposed to manual creation.  So, if there is another language with this capability, please suggest it (I did see the <a href="https://stackoverflow.com/questions/9071084/polar-contour-plot-in-matplotlib-best-modern-way-to-do-it">Python example</a>).</p>
  
  <h3>EDIT</h3>
  
  <p>Regarding the suggestion using ggplot2 - I can't seem to get the geom_tile routine to plot interpolated data in polar_coordinates.  I have included code below that illustrates where I'm at.  I can plot the original in Cartesian and polar, but I can only get the interpolated data to plot in Cartesian.  I can plot the interpolation points in polar using geom_point, but I can't extend that approach to geom_tile.  My only guess was that this is related to data order - i.e. geom_tile is expecting sorted/ordered data - I've tried every  iteration I can think of sorting the data into ascending/descending azimuth and zenith with no change.</p>
  
  <pre><code>## Libs
  library(akima)
  library(ggplot2)
  
  ## Sample data in az/el(zenith)
  tmp = seq(5,355,by=10)
  geoms &lt;- 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&lt;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&lt;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()
  </code></pre>
  
  <h3>Final Results</h3>
  
  <p>I finally took code from the accepted answer (base graphics) and updated the code.  I added a thin plate spline interpolation method, an option to extrapolate or not, data point overlays, and the ability to do continuous colors or segmented colors for the interpolated surface.  See the examples below.</p>
  
  <pre><code>PolarImageInterpolate &lt;- 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 &lt;- seq(-outer.radius, outer.radius, length.out = spatial_res)
  # interpolate the data
      if (interp.type ==1 ){
      Interp &lt;- akima:::interp(x = x, y = y, z = z, 
                      extrap = extrapolate, 
                      xo = minitics, 
                      yo = minitics, 
                      linear = FALSE)
      Mat &lt;- 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 &lt;- matrix(minitics, ncol = spatial_res, nrow = spatial_res) 
  Mat[!sqrt(markNA ^ 2 + t(markNA) ^ 2) &lt; outer.radius] &lt;- NA 
  
      ### Set contour_breaks based on requested source
      if ((length(contour_breaks_source == 1)) &amp; (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)) &amp; (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) &amp; (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) &amp; (col_breaks_source[1] == 1))
          {zlim=c(min(z,na.rm=TRUE),max(z,na.rm=TRUE))}
      else if ((length(col_breaks_source) == 1) &amp; (col_breaks_source[1] == 2))
          {zlim=c(min(Mat,na.rm=TRUE),max(Mat,na.rm=TRUE))}
      else if ((length(col_breaks_source) == 2) &amp; (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&lt;zlim[1])]=zlim[1]
      Mat_plot[which(Mat_plot&gt;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 &lt;- contourLines(x = minitics, y = minitics, Mat, levels = contour_breaks)
      A &lt;- 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 &lt;- function(radians){
          matrix(c(cos(radians), sin(radians), -sin(radians), cos(radians)), ncol = 2)
      }    
  
      circle &lt;- function(x, y, rad = 1, nvert = 500){
          rads &lt;- seq(0,2*pi,length.out = nvert)
          xcoords &lt;- cos(rads) * rad + x
          ycoords &lt;- sin(rads) * rad + y
          cbind(xcoords, ycoords)
      }
  
      # draw circles
      if (missing(circle.rads)){
          circle.rads &lt;- pretty(c(0,outer.radius))
      }
  
      for (i in circle.rads){
          lines(circle(0, 0, i), col = "#66666650")
      }
  
      # put on radial spoke axes:
      axis.rads &lt;- c(0, pi / 6, pi / 3, pi / 2, 2 * pi / 3, 5 * pi / 6)
      r.labs &lt;- c(90, 60, 30, 0, 330, 300)
      l.labs &lt;- c(270, 240, 210, 180, 150, 120)
  
      for (i in 1:length(axis.rads)){ 
          endpoints &lt;- 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 &lt;- c(RMat(axis.rads[i]) %*% matrix(c(1.1, 0, -1.1, 0) * outer.radius, ncol = 2))
          lab1 &lt;- bquote(.(r.labs[i]) * degree)
          lab2 &lt;- 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 &lt;- seq(-outer.radius, outer.radius, length = contour_levels+ 1)
      # #ylevs &lt;- 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)
      }
  }
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/XawPX.png" alt="enter image description here">
  <img src="https://i.stack.imgur.com/Fj8wW.png" alt="enter image description here">
  <img src="https://i.stack.imgur.com/urbh0.png" alt="enter image description here"></p>