Develop a modified version of stat_contour

Question

I'm ultimately trying to plot contour plots, or "raster plots", of irregular datasets - a rather common question of course. Many solutions propose to interpolate the data first, and then plot it, for instance here : Plotting contours on an irregular grid amongst other - or in fact, the man page at https://ggplot2.tidyverse.org/reference/geom_contour.html

However, for convenience I'm trying to wrap it into a new stat.

I managed to get something that works for geom_raster, simply lifting the interpolation code from the example in the manual:

require(akima)
StatInterpRaster <- ggproto("StatInterpRaster", Stat,
                           compute_group = function(data, scales) {
                             
                             ii<-akima::interp(x = data$x,
                                               y = data$y,
                                               z = data$fill)
                             
                             data.out <- tibble(x = rep(ii$x, nrow(ii$z)),
                                           y = rep(ii$y, each = ncol(ii$z)),
                                           fill = as.numeric(ii$z) )
                         
                             return(data.out)
                           },
                           
                           required_aes = c("x", "y", "fill")
)

stat_interp_raster<- function(mapping = NULL, data = NULL, geom = "contour",
                              position = "identity", na.rm = FALSE, show.legend = NA,
                              inherit.aes = TRUE, ...) {
  layer(
    stat = StatInterpRaster, data = data, mapping = mapping, geom = geom,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}

which works as expected:

ee <- tibble (x=rnorm(50),y=rnorm(50),z=x*y)
ee %>% ggplot() + geom_raster(aes(x=x,y=y,fill=z),stat=StatInterpRaster)

I would now trying to achieve the same thing with contours. Naively I tried

StatInterpContour <- ggproto("StatInterpContour", Stat,
                            compute_group = function(data, scales) {
                              
                              ii<-akima::interp(x = data$x,
                                                y = data$y,
                                                z = data$z)
                              
                              data.out <- tibble(x = rep(ii$x, nrow(ii$z)),
                                                 y = rep(ii$y, each = ncol(ii$z)),
                                                 z = as.numeric(ii$z) )
                              
                              #StatContour(data.out)
                              return(data.out)
                            },
                            
                            required_aes = c("x", "y", "z")
)

stat_interp_contour<- function(mapping = NULL, data = NULL, geom = "contour",
                              position = "identity", na.rm = FALSE, show.legend = NA,
                              inherit.aes = TRUE, ...) {
  layer(
    stat = StatInterpContour, data = data, mapping = mapping, geom = geom,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, ...)
  )
}

which is essentially the same as above. However it does not produce the expected result :

 ee %>% ggplot() + geom_contour(aes(x=x,y=y,z=z),stat=StatInterpContour)

In retrospect, this is not surprising. My stat is generating a regular data array, with neatly ordered values in x and y, but nowhere am I generating the actual lines. The contour lines are more complicated, seem to be generated by xyz_to_isolines in stat_contour (cf. https://github.com/tidyverse/ggplot2/blob/main/R/stat-contour.r , line 97 as of today).

I could copy the relevant code in stat-contour.r, but it seems to me that it is a waste of effort and it would be better to simply pass my result to stat_contour, that already does the job: it generates contour lines from an object of that shape. So apparently I "just" have to call StatContour (or friends) somewhere in my StatInterpContour function compute_group -- but how ?

Thanks !

Answer 1

You are right that you shouldn't need to copy code over from StatContour. Instead, make your ggproto class inherit from StatContour. Prepare the data then pass it, along with all necessary parameters, to the compute_group function from StatContour

StatInterpContour <- ggproto("StatInterpRaster", StatContour,
  compute_group = function(data, scales, z.range, bins = NULL, binwidth = NULL, 
    breaks = NULL, na.rm = FALSE)  {
                             
    ii<-akima::interp(x = data$x,
                      y = data$y,
                      z = data$z)
    data <-  tibble(x = rep(ii$x, nrow(ii$z)),
             y = rep(ii$y, each = ncol(ii$z)),
             z = as.numeric(ii$z), group = 1) 
    StatContour$compute_group(data, scales, z.range, 
                  bins, binwidth, breaks, na.rm)
    
    },
    required_aes = c("x", "y", "z")
)

This requires a little modification of your user-facing function:

stat_interp_contour<- function(mapping = NULL, data = NULL, geom = "contour",
                              position = "identity", na.rm = FALSE, 
                              show.legend = NA,
                              inherit.aes = TRUE, bins = NULL, binwidth = NULL, 
    breaks = NULL,  ...) {
  layer(
    stat = StatInterpContour, data = data, mapping = mapping, geom = geom,
    position = position, show.legend = show.legend, inherit.aes = inherit.aes,
    params = list(na.rm = na.rm, bins = bins, binwidth = binwidth, 
    breaks = breaks,  ...)
  )
}

But should now work without as expected. Here, I've plotted it along with the original points coloured according to their z value to show that the contours try to approximate the level of the points:

ee %>% 
  ggplot(aes(x, y)) + 
  geom_point(aes(color = z), size = 3) +
  stat_interp_contour(aes(z = z, color = after_stat(level))) +
  scale_color_viridis_c()