Producing a raster from a ggplot2 plot

Question

I am working with a plot that I need to somehow apply a unique regression to. The plot is a series of points placed on log-log axes. For anyone familiar with hydrology, I'm trying to graphically produce a transient Theis solution, so my regression needs to be in the form of an exponential integral (see this Wikipedia page).

It's important to note here that the exponential integral has its own set of axes values that are independent of the initial plot. This is how a solution is pulled from the data, but it introduces a number of problems when trying to reproduce it in R.

I've come up with a few ideas to deal with this (beyond using a pencil and paper), but have hit some small snags with each of the approaches. I'd appreciate any insight into workarounds for any of these solutions:

1. Adding a regression to the plot using stat_smooth, arbitrarily selecting a point on the regression based on its slope, then correlating that to the same slope of a new plot that has the appropriate axes for the exponential integral. The issue here is that I am not sure how to use stat_smooth with a formula other than variations on y ~ x, y ~ poly(x, 2), and y ~ log(x). The regression would call for a combination of the poly and log functions, but R hiccups when I try to do this. For example:

   stat_smooth(data = example, method = "lm",
       formula = y ~ log(x) + x - poly(x, 2)/4 + poly(x, 3)/18 - ...
   

2. I've also tried to plot the data and the exponential integral as two separate plots, then overlay the two plots based on where I arbitrarily believe the best-fit to be (this may be the easier method, and is plenty accurate enough for my purposes). To facilitate the overlay, I've stripped the exponential integral plot of its background, axes, gridlines, etc. and left only a horizontal and vertical line with labels to indicate a point with certain values on the curve. If I can place this on top of the other plot (assuming, of course, that both plots maintain the same size ratio), I figure I can nudge the regression around until it lines up where I think it should, then read off its corresponding values based on the presence of the horizontal and vertical labeled lines.

   I've read a little bit about annotation_raster, and think it may be an appropriate way for me to treat the regression as an overlayed image. My issue, though, is with converting a ggplot plot to a raster first. as.raster() produces the following error:

   Error in as.raster(raster) : 
     error in evaluating the argument 'x' in selecting a method for function 'as.raster':
   Error in UseMethod("as.raster") : 
     no applicable method for 'as.raster' applied to an object of class "c('gg', 'ggplot')"
   

   A similar error appears if I try to use the raster package and convert the plot using the raster() function. Is there a simple way to do this?

Here's some reproducible code:

library(ggplot2)

Data = data.frame(matrix(
    # Elapsed_sec  Drawdown_ft
    c(20,  0.0038,
      40,  0.0094,
      60,  0.017,
      80,  0.0283,
      100, 0.0358,
      120, 0.0415,
      140, 0.049,
      160, 0.0528,
      180, 0.0548,
      200, 0.0567), nrow = 10, ncol = 2, byrow = TRUE))
colnames(Data) = c("Elapsed_sec", "Drawdown_ft")

Integral = data.frame(matrix(
    # u     W_u
    c(1e-3, 6.33,
      5e-3, 4.73,
      1e-2, 4.04,
      5e-2, 2.47,
      1e-1, 1.82,
      5e-1, 0.56,
      1e0,  0.219,
      2e0,  0.049,
      3e0,  0.013,
      4e0,  0.0038,
      5e0,  0.0011,
      6e0,  0.00036), nrow = 12, ncol = 2, byrow = TRUE))
colnames(Integral) = c("u", "W_u")

# Plot exponential integral (Theis curve)
Tcurve = ggplot(Integral, aes(1/u, W_u)) + geom_line() +
    scale_x_log10(limits = c(10^-1, 10^3), breaks = c(10^-1, 10^0,  10^1,  10^2, 10^3)) +
    scale_y_log10(limits = c(10^-3, 10^1), breaks = c(10^-3, 10^-2, 10^-1, 10^0, 10^1)) +
    xlab("1/u") + ylab("W(u)") + coord_equal() +
    geom_hline(aes(yintercept = 0.219), linetype = 2) +
    geom_vline(aes(xintercept = 1),     linetype = 2) +
    geom_text(color = "black", size = 3, aes(x = 0.3, y = 0.3,  label = "W(u) = 0.219")) +
    geom_text(color = "black", size = 3, aes(x = 0.8, y = 0.01, label = "u = 1"), angle = 90) +
    theme(line = element_blank(), text = element_blank(), panel.background = element_rect(fill = NA))

# Plot drawdown data
plot = ggplot(Data, aes(Elapsed_sec, Drawdown_ft)) + geom_point(alpha = 0.5, size = 1) +
    scale_x_log10(limits = c(10^0,  10^4), breaks = c(10^0,  10^1,  10^2,  10^3, 10^4)) +
    scale_y_log10(limits = c(10^-3, 10^1), breaks = c(10^-3, 10^-2, 10^-1, 10^0, 10^1)) +
    xlab("Elapsed Time (sec)") + ylab("Drawdown (ft)") + coord_equal() + theme_bw()

I'd upload images, but I currently lack the minimum reputation. The first plot shows the exponential integral, while the second shows the data that I am trying to fit the integral to.

Answer 1

It has been five years since this question posted, but I thought I would share my thoughts. You may convert a ggplot to a raster object as follows:

1. Save the ggplot plot into an image file of any format (I would use tiff) using ggsave
2. Create a raster object for the saved image using raster (or stack for colored images)

The implementation of the above is as follows:

# You may need to remove the margins from the plot (i.e., keep the earth raster and the routes only)
plot <- plot+ 
  theme(    
        axis.ticks=element_blank(), 
        axis.text.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.title.x=element_blank(), 
        axis.title.y=element_blank(),
        plot.margin = unit(c(0, 0, 0, 0), "null"),
        legend.position = 'none'
       ) +
       labs(x=NULL, y=NULL)

# Save the plot
ggsave(plot=plot, "my_ggplot.tiff", device = "tiff")

# Create a StackedRaster object from the saved plot
raster <- raster("my_ggplot.tiff") # OR stack("my_ggplot.tiff") for colored images

# Get the GeoSpatial Components
lat_long <- ggplot_build(plot)$layout$panel_params[[1]][c("x.range","y.range")] 

# Supply GeoSpatial  data to the StackedRaster 
extent(raster) <- c(lat_long$x.range,lat_long$y.range)
projection(raster) <- CRS("+proj=longlat +datum=WGS84")


# Then, do whatever you want with the raster object here

Hope it helps.