Obtain Function from Akima::Interp() Matrix

Obtain function from akima::interp() matrix

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.

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)

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()

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)

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()

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 ?

Interp() Function Akima Plots Confidence/Error Terms

The akima interpolation model is not a statistical model, its an interpolation, so there's no model for the variance which is what you'd need to get any measure of "confidence".

Interpolation (akima) omits part of data when x/y contains duplicate elements

It turns out that the error arose from plot_ly(). Apparently, the z-matrix cannot be passed straight through from interp() to plot_ly(), as the axis become erroneously passed through to the graph. Hence, the interpolated z-matrix needs to be transformed.

If you use these two functions in combination, ensure to carry out the transformation of z as shown below:

mat = interp(x,y,z, duplicate = 'mean')

x = mat$x
y = mat$y
z = matrix(mat$z, nrow = length(mat$y), byrow = TRUE)

plot_ly(x, y ,z, type = 'surface')

How to improve interp with akima?

From the documentation for interp, regarding the xo argument:

vector of x-coordinates of output grid. The default is 40 points
evenly spaced over the range of x.

So the x in fld may not include the x-value corresponding to your maximum or minimum a1 values.

By increasing the number of points for xo and yo, the max z value interpolated will get closer to the actual max. For example,

set.seed(100)
x <- rnorm(1398)
y <- rnorm(1398)
a1 <- rnorm(1398)
data <- data.frame(x, y, a1)
fld <- with(data, interp(x,y,a1))
fld2 <- with(data, interp(x,y,a1, xo=seq(min(x), max(x), length=1000), yo=seq(min(y), max(y), length=1000)))

max(a1, na.rm=TRUE)
#[1] 2.949
max(fld$z, na.rm=TRUE)
#[1] 2.481
max(fld2$z, na.rm=TRUE)
#[1] 2.902

Furthermore, if you want the interpolated z to include your max and min a1, add corresponding x and y values to xo and yo. For example, this is how you would get it to include the max value of a1.

max.a1.x <- x[which.max(a1)]
max.a1.y <- y[which.max(a1)]
# these have to be sorted, since filled.contour will expect them to be.
xo <- sort(c(seq(min(x), max(x), length=40), max.a1.x))
yo <- sort(c(seq(min(y), max(y), length=40), max.a1.y))

fld3 <- with(data, interp(x,y,a1, xo=xo, yo=yo))
filled.contour(x=fld3$x,y=fld3$y,z=fld3$z,
  color.palette=colorRampPalette(c("white", "blue")))

max(fld3$z, na.rm=TRUE)
# [1] 2.949

---

Related Topics