将不规则网格插值到规则网格答案

【问题标题】：Interpolate irregular grid to regular grid将不规则网格插值到规则网格
【发布时间】：2014-03-08 14:41:21
【问题描述】：

我有一个不规则网格，我需要将其转换为规则网格，以利用 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 年。所以我想我的新问题只是整理上面的代码来执行聚合方法。

【问题讨论】：

真正的数据源是什么？ netcdf中的东西？值得研究
@mdsumner 时间与年龄，两个维度都有年的单位。 z 是随年龄变化的概率。

标签： r interpolation

【解决方案1】：

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

我使用多核技术解决了这个问题，所以如果您有一个多核处理器，请考虑类似于以下代码 sn-p 的内容：

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]] ])) )} )

这是从我为在图像上执行旋转“漩涡”而编写的一个函数中提取的，所以那里有很多你不需要的东西。

【讨论】：