从 akima::interp() 矩阵获取函数
Obtain function from akima::interp() matrix
使用 interp 函数（Akima 包），可以绘制对应于数据集的双变量插值的曲面，请参见下面的示例（来自 interp 文档）：
library(rgl)
data(akima)
# data visualisation
rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red")
rgl.bbox()
# bivariate linear interpolation
# interp:
akima.li <- interp(akima$x, akima$y, akima$z,
xo=seq(min(akima$x), max(akima$x), length = 100),
yo=seq(min(akima$y), max(akima$y), length = 100))
# interp surface:
rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5))
但是，输出只是描述一组点的列表，而不是一般函数。
问题：有什么方法可以得到一个函数z = f(x,y)来匹配之前得到的曲面吗？我知道它可以使用 interp(akima$x, akima$y, akima$z, xo=A, yo=B)，但速度很慢。
在二维中，approxfun() 函数可以完成这项工作，但我找不到多参数插值的等价物。
如果你想要一个线性插值使得表面穿过所有点，你将无法使用函数进行插值z = f(x,y)，除非数据集已经通过这种函数进行了模拟。
如果您正在寻找与您的点集相匹配的函数 z=f(x,y)，则必须使用 GLM 或 GAM 等构建模型。但是，这会导致表面不会跨越所有点数据，并且会存在一些残差。
因为我过去常常使用空间数据集，这意味着 x 和 y 坐标与 z 观测值，我将以这种方式为您提供一些线索。
首先，我准备了一个插值数据集：
library(rgl)
library(akima)
library(dplyr)
library(tidyr)
data(akima)
data.akima <- as.data.frame(akima)
# data visualisation
rgl.spheres(akima$x, akima$z , akima$y,0.5,color="red")
rgl.bbox()
# Dataset for interpolation
seq_x <- seq(min(akima$x) - 1, max(akima$x) + 1, length.out = 20)
seq_y <- seq(min(akima$y) - 1, max(akima$y) + 1, length.out = 20)
data.pred <- dplyr::full_join(data.frame(x = seq_x, by = 1),
data.frame(y = seq_y, by = 1)) %>%
dplyr::select(-by)
然后，我用你的akima插值函数：
# bivariate linear interpolation
# interp:
akima.li <- interp(akima$x, akima$y, akima$z,
xo=seq_x,
yo=seq_y)
# interp surface:
rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5))
rgl.spheres(akima$x, akima$z , akima$y,0.5,color="red")
rgl.bbox()
使用栅格
从现在开始，如果您想获取某些特定点的插值信息，您可以重新使用 interp 函数或决定使用光栅化图像。使用栅格，您可以提高分辨率，并获得任何空间位置信息数据。
# Using rasters
library(raster)
r.pred <- raster(akima.li$z, xmn = min(seq_x), xmx = max(seq_x),
ymn = min(seq_y), ymx = max(seq_y))
plot(r.pred)
## Further bilinear interpolations
## Double raster resolution
r.pred.2 <- disaggregate(r.pred, fact = 2, method = "bilinear")
plot(r.pred.2)
空间插值（反距离插值或克里金法）
在考虑空间插值时，我首先想到的是克里金法。这将使您的表面平滑，因此它不会跨越每个数据点。
# Spatial inverse distance interpolation
library(sp)
library(gstat)
# Transform data as spatial objects
data.akima.sp <- data.akima
coordinates(data.akima.sp) <- ~x+y
data.pred.sp <- data.pred
coordinates(data.pred.sp) <- ~x+y
# Inverse distance interpolation
# idp is set to 2 as weight for interpolation is :
# w = 1/dist^idp
# nmax is set to 3, so that only the 3 closest points are used for interpolation
pred.idw <- idw(
formula = as.formula("z~1"),
locations = data.akima.sp,
newdata = data.pred.sp,
idp = 2,
nmax = 3)
data.spread.idw <- data.pred %>%
select(-pred) %>%
mutate(idw = pred.idw$var1.pred) %>%
tidyr::spread(key = y, value = idw) %>%
dplyr::select(-x)
surface3d(seq_x, seq_y, as.matrix(data.spread.idw), col = "green")
rgl.spheres(akima$x, akima$y , akima$z, 0.5, color = "red")
rgl.bbox()
使用 gam 或 glm 进行插值
但是，如果你想找到像z = f(x,y)这样的公式，你应该使用自由度高的GLM或GAM，这取决于你希望看到的光滑度。另一个优点是您可以添加其他协变量，而不仅仅是 x 和 y。该模型需要安装 x/y 交互。
这里有一个简单的 GAM 平滑示例：
# Approximation with a gam model
library(mgcv)
gam1 <- gam(z ~ te(x, y), data = data.akima)
summary(gam1)
plot(gam1)
data.pred$pred <- predict(gam1, data.pred)
data.spread <- tidyr::spread(data.pred, key = y, value = pred) %>%
dplyr::select(-x)
surface3d(seq_x, seq_y, as.matrix(data.spread), col = "blue")
rgl.spheres(akima$x, akima$y , akima$z, 0.5, color = "red")
rgl.bbox()
这个答案对您来说方向正确吗？
使用 interp 函数（Akima 包），可以绘制对应于数据集的双变量插值的曲面，请参见下面的示例（来自 interp 文档）：
library(rgl)
data(akima)
# data visualisation
rgl.spheres(akima$x,akima$z , akima$y,0.5,color="red")
rgl.bbox()
# bivariate linear interpolation
# interp:
akima.li <- interp(akima$x, akima$y, akima$z,
xo=seq(min(akima$x), max(akima$x), length = 100),
yo=seq(min(akima$y), max(akima$y), length = 100))
# interp surface:
rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5))
但是，输出只是描述一组点的列表，而不是一般函数。
问题：有什么方法可以得到一个函数z = f(x,y)来匹配之前得到的曲面吗？我知道它可以使用 interp(akima$x, akima$y, akima$z, xo=A, yo=B)，但速度很慢。
在二维中，approxfun() 函数可以完成这项工作，但我找不到多参数插值的等价物。
如果你想要一个线性插值使得表面穿过所有点，你将无法使用函数进行插值z = f(x,y)，除非数据集已经通过这种函数进行了模拟。
如果您正在寻找与您的点集相匹配的函数 z=f(x,y)，则必须使用 GLM 或 GAM 等构建模型。但是，这会导致表面不会跨越所有点数据，并且会存在一些残差。
因为我过去常常使用空间数据集，这意味着 x 和 y 坐标与 z 观测值，我将以这种方式为您提供一些线索。
首先，我准备了一个插值数据集：
library(rgl)
library(akima)
library(dplyr)
library(tidyr)
data(akima)
data.akima <- as.data.frame(akima)
# data visualisation
rgl.spheres(akima$x, akima$z , akima$y,0.5,color="red")
rgl.bbox()
# Dataset for interpolation
seq_x <- seq(min(akima$x) - 1, max(akima$x) + 1, length.out = 20)
seq_y <- seq(min(akima$y) - 1, max(akima$y) + 1, length.out = 20)
data.pred <- dplyr::full_join(data.frame(x = seq_x, by = 1),
data.frame(y = seq_y, by = 1)) %>%
dplyr::select(-by)
然后，我用你的akima插值函数：
# bivariate linear interpolation
# interp:
akima.li <- interp(akima$x, akima$y, akima$z,
xo=seq_x,
yo=seq_y)
# interp surface:
rgl.surface(akima.li$x,akima.li$y,akima.li$z,color="green",alpha=c(0.5))
rgl.spheres(akima$x, akima$z , akima$y,0.5,color="red")
rgl.bbox()
使用栅格
从现在开始，如果您想获取某些特定点的插值信息，您可以重新使用 interp 函数或决定使用光栅化图像。使用栅格，您可以提高分辨率，并获得任何空间位置信息数据。
# Using rasters
library(raster)
r.pred <- raster(akima.li$z, xmn = min(seq_x), xmx = max(seq_x),
ymn = min(seq_y), ymx = max(seq_y))
plot(r.pred)
## Further bilinear interpolations
## Double raster resolution
r.pred.2 <- disaggregate(r.pred, fact = 2, method = "bilinear")
plot(r.pred.2)
空间插值（反距离插值或克里金法）
在考虑空间插值时，我首先想到的是克里金法。这将使您的表面平滑，因此它不会跨越每个数据点。
# Spatial inverse distance interpolation
library(sp)
library(gstat)
# Transform data as spatial objects
data.akima.sp <- data.akima
coordinates(data.akima.sp) <- ~x+y
data.pred.sp <- data.pred
coordinates(data.pred.sp) <- ~x+y
# Inverse distance interpolation
# idp is set to 2 as weight for interpolation is :
# w = 1/dist^idp
# nmax is set to 3, so that only the 3 closest points are used for interpolation
pred.idw <- idw(
formula = as.formula("z~1"),
locations = data.akima.sp,
newdata = data.pred.sp,
idp = 2,
nmax = 3)
data.spread.idw <- data.pred %>%
select(-pred) %>%
mutate(idw = pred.idw$var1.pred) %>%
tidyr::spread(key = y, value = idw) %>%
dplyr::select(-x)
surface3d(seq_x, seq_y, as.matrix(data.spread.idw), col = "green")
rgl.spheres(akima$x, akima$y , akima$z, 0.5, color = "red")
rgl.bbox()
使用 gam 或 glm 进行插值
但是，如果你想找到像z = f(x,y)这样的公式，你应该使用自由度高的GLM或GAM，这取决于你希望看到的光滑度。另一个优点是您可以添加其他协变量，而不仅仅是 x 和 y。该模型需要安装 x/y 交互。
这里有一个简单的 GAM 平滑示例：
# Approximation with a gam model
library(mgcv)
gam1 <- gam(z ~ te(x, y), data = data.akima)
summary(gam1)
plot(gam1)
data.pred$pred <- predict(gam1, data.pred)
data.spread <- tidyr::spread(data.pred, key = y, value = pred) %>%
dplyr::select(-x)
surface3d(seq_x, seq_y, as.matrix(data.spread), col = "blue")
rgl.spheres(akima$x, akima$y , akima$z, 0.5, color = "red")
rgl.bbox()
这个答案对您来说方向正确吗？