从ggplot2绘图生成栅格图像

Question

我正在处理一个绘图问题，需要对其应用独特的回归分析。这个图是一系列点在对数对数坐标轴上的展示。对于熟悉水文学的人来说，我试图以图像方式产生瞬态Theis解，因此我的回归必须采用指数积分的形式（请参见维基百科网页）。

需要注意的是，指数积分具有自己独立于初始图的一组坐标轴值。这是从数据中提取解决方案的方法，但在尝试在R中重现它时会引入许多问题。

我已经想出了一些应对措施（除了使用铅笔和纸），但每种方法都遇到了一些小问题。如果您能提供任何解决方案的见解，我将非常感激。

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?

这里是一些可复现的代码：

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

我想上传图片，但是目前我的声望还不够。第一个图显示了指数积分，而第二个图显示了我试图将该积分拟合到的数据。

- hfisch

---

你能否提供一个可重现的例子（即一些虚拟数据和可重现的代码）。你不能从ggplot对象强制转换为raster - 毫无疑问，会有一个简单的解决方案，但如果没有可重现的例子，建议一个将非常困难。 - mnel

@mnel 已添加可重现的代码。我已大大减少了数据集，如果添加更多数据点有帮助，请告诉我。 - hfisch

Answer 1

自从这个问题发布已经五年了，但我想分享一下我的想法。您可以按照以下步骤将ggplot转换为raster对象：

1. 使用ggsave将ggplot图形保存为任何格式的图像文件（我会使用tiff）
2. 使用raster（或对于彩色图像，使用stack）为保存的图像创建raster对象

实现如下：

# 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

希望能对您有所帮助。

- Taher A. Ghaleb