从akima :: interp()矩阵获取函数 [英] Obtain function from akima::interp() matrix
问题描述
使用interp函数(Akima软件包),可以绘制与数据集的双变量插值相对应的曲面,请参见下面的示例(来自interp文档):
Using the interp function (Akima package), it is possible to draw the surface corresponding to the bivariate interpolation of a data set, see example below (from interp documentation):
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))
但是,输出仅是描述一组点的列表,而不是一般函数。
However, the output is only a list describing a set of points, not a general function.
问题:是否可以通过任何方法获得与先前函数匹配的函数z = f(x,y)获得的表面?我知道它可以使用interp(akima $ x,akima $ y,akima $ z,xo = A,yo = B),但是它非常慢。
Question: is there any method to obtain a function z = f(x,y) that matches the previously obtained surface ? I know that it works using interp(akima$x, akima$y, akima$z, xo=A, yo=B), but it is very slow.
在二维中,roximfun()函数可以完成此工作,但是我找不到等效的多个参数插值。
In two dimensions, the approxfun() function would do the job, but I could not find the equivalent for multiple parameters interpolation.
推荐答案
如果要进行线性插值以使曲面与所有点交叉,则无法使用函数 z = f(x,y)
进行插值,除非数据集已通过这种函数进行了模拟。
如果您正在寻找与点集匹配的函数 z = f(x,y)
,您将必须使用例如GLM或GAM建立模型。但是,这导致曲面不会与所有点数据交叉,并且会存在一些残差。
If you want a linear interpolation so that the surface cross all points, you will not be able to interpolate with a function z = f(x,y)
, except if the dataset has been simulated through this kind of function.
If you are looking for a function z=f(x,y)
that matches your point set, you will have to build a model with GLM or GAM for instance. However, this induces that the surface will not cross all points data and there will be some residuals.
在我过去使用空间数据集时,这意味着x和y坐标与z观测值一致,我将以这种方式为您提供一些线索。
As I use to work with spatial datasets, which means x and y coordinates with a z observation, I will give you some clues in this way.
首先,我准备一个用于插值的数据集:
First, I prepare a dataset for interpolation:
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插值函数:
Then, I use your akima interpolation function:
# 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()
< h1>使用栅格
从现在开始,如果要获取某些特定点的插值信息,可以重新使用 interp
功能或决定使用栅格化图像。使用栅格,您便可以提高分辨率,并获取任何空间位置信息数据。
Using rasters
From now, if you want to get interpolated information on some specific points, you can re-use interp
function or decide to work with a rasterized image. Using rasters, you are then able to increase resolution, and get any spatial position information data.
# 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)
空间插值(逆距离插值或kriging)
在考虑空间插值时,我首先想到了kriging。这样会平滑您的表面,因此不会跨越每个数据点。
Spatial interpolation (inverse distance interpolation or kriging)
When thinking in spatial for interpolation, I first think about kriging. This will smooth your surface, thus it will not cross every data points.
# 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。该模型需要适应ax / y交互。
下面是一个具有简单GAM平滑的示例:
Interpolate using gam or glm
However, if you want to find a formula like z = f(x,y)
, you should use GLM or GAM with high degrees of freedom depending on the smooth you hope to see. Another advantage is that you can add other covariates, not only x and y. The model needs to be fitted with a x/y interaction.
Here an example with a simple GAM smooth:
# 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()
这个答案对您来说是正确的方向吗?
Does this answer goes in the right direction for you ?
这篇关于从akima :: interp()矩阵获取函数的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!