从akima :: interp（）矩阵获取函数

本文介绍了从akima :: interp（）矩阵获取函数的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述 使用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 visualisationrgl.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 visualisationrgl.spheres(akima$x, akima$z , akima$y,0.5,color="red")rgl.bbox()# Dataset for interpolationseq_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() 使用栅格 从现在开始，如果要获取某些特定点的插值信息，可以重新使用 interp 功能或决定使用栅格化图像。使用栅格，您便可以提高分辨率，并获取任何空间位置信息数据。Using rastersFrom 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 resolutionr.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 interpolationlibrary(sp)library(gstat)# Transform data as spatial objectsdata.akima.sp <- data.akimacoordinates(data.akima.sp) <- ~x+ydata.pred.sp <- data.predcoordinates(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 interpolationpred.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 glmHowever, 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 modellibrary(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（）矩阵获取函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！