r - 将不规则网格插入到规则网格

Question

我有一个不规则的网格，我需要将其转换为规则网格才能利用 image的 useRaster=TRUE图形设备的选项。我可以通过将不规则网格转换为点来小规模地执行此操作，然后使用 akima 的 interp 对这些点进行插值。 .然而，这会随着尺寸的增加而可怕地扩展，所以我正在寻找选择。

首先，这里是小比例 (5x10) 示例，其中只有 x 维度是不规则的:

nx <- 5
ny <- 10
si <- list()  # irregular surface
si$x <- cumsum(runif(nx) * 10) + 100
si$y <- seq(20, 50, length.out=ny)
si$z <- matrix(rnorm(nx * ny), ncol=ny)
image(si)

和双线性插值结果:

sr_x <- seq(min(si$x), max(si$x), length.out=nx * 5)
sr_y <- si$y  # this dimension is already regular
require(akima)  # interpolate from points repeated off irregular grid
sr <- interp(rep(si$x, length(si$y)), rep(si$y, each=length(si$x)), si$z,
             xo=sr_x, yo=sr_y)
image(sr, useRaster=TRUE)

但是，如果使用更大尺寸的不规则网格(例如 nx <- 50; ny <- 100 )，过程会很慢。是否有可以加快该过程的库或技术？

---

更新和可能的解决方案。数据描述了时间与时间(均以年为单位)，其中不规则维度的时间步长介于 0.5 天到 30 天之间，第二个时间轴具有等间距365天间隔。由于沿不规则轴的间距要小得多，因此插值将不起作用。因此，平滑或聚合方法会产生更好的结果。

更真实的数据场景，展现更精细的不规则维度:

nx <- 200
ny <- 10
si <- list()  # irregular surface
si$x <- cumsum(runif(nx, 0.5, 30) / 365)
si$y <- 1:ny
si$z <- matrix(rnorm(nx * ny), ncol=ny)
image(si)

还有一些非常粗略的聚合意味着:

dx <- 1/12  # 1 month spacing along x-axis
sr <- list()  # regular surface
sr$x <- seq(min(si$x), max(si$y), dx)  # equal-width breaks
nsrx <- length(sr$x)
sr$y <- si$y  # this dimension is already regular
sr$z <- matrix(nrow=length(sr$x), ncol=length(sr$y))
# Classify irregular dimension
si_xc <- cut(si$x, sr$x, include.lowest=TRUE, labels=FALSE)
# Aggregate means from irregular to regular dimension
for(xi in seq_len(nsrx))
    sr$z[xi,] <- apply(si$z[si_xc == xi, , drop=FALSE], 2, mean)
image(sr, zlim=range(si$z), useRaster=TRUE)

这似乎可以解决问题，并且可以在每个维度上具有 100 年的更大的数据集上进行扩展。所以我想我的新问题只是整理上面的代码以执行聚合方法。

Answer 1

有几个带有“kriging”工具的软件包，基本上就是您想要的。但是，我不知道它是否会比 akima::interp 做的更快。

我使用多核技术解决了这个问题，所以如果您有一个多核处理器，请考虑类似于以下代码片段的东西:

picbits <- clusterApply( myclus, 1:length(picsec) , function(j) { gc(); 
akima::interp(newx[picsec[[j]] ], newy[picsec[[j]] ], picture[picsec[[j]] ], 

xo=trunc(min(newx[picsec[[j]] ])):trunc(max(newx[picsec[[j]] ])), 

yo=trunc(min(newy[picsec[[j]] ])):trunc(max(newy[picsec[[j]] ])) )} )

这是从我编写的用于在图像上执行旋转“漩涡”的函数中提取的，因此那里有很多你不需要的东西。