Search code examples
rinterpolationbilinear-interpolation

Error when doing bilinear interpolation with `interp2 {pracma}`; any better way for 2D interpolation?


I am trying to perform 2d interpolation on a table named vol_coarse.

install.packages("install.load")
install.load::load_package("pracma", "data.table")

vol_coarse <- data.table(V1 = c(3 / 8, 1 / 2, 3 / 4, 1, 1 + 1 / 2, 2, 3, 6),
V2 = c(0.50, 0.59, 0.66, 0.71, 0.75, 0.78, 0.82, 0.87),
V3 = c(0.48, 0.57, 0.64, 0.69, 0.73, 0.76, 0.80, 0.85),
V4 = c(0.44, 0.53, 0.60, 0.65, 0.69, 0.72, 0.76, 0.81))
setnames(vol_coarse, c("Maximum size of aggregate (in)", "2.40", "2.60", "2.80"))

x <- vol_coarse[, 2][[1]]

y <- as.numeric(colnames(vol_coarse[, 2:ncol(vol_coarse)]))

z <- meshgrid(x, y)

xp <- 3 / 4

yp <- 2.70

interp2(x = x, y = y, Z = z, xp = xp, yp = yp, method = "linear")

This is the error message that is returned:

Error: is.numeric(Z) is not TRUE

I read in ?interp2 that:

length(x) = nrow(Z) = 8 and length(y) = ncol(Z) = 3 must be satisfied.

How can I create a matrix that is 8 by 3 so that I can use interp2?

Or is there a better way to perform this type of interpolation?

Thank you.


Solution

  • If I don't get you wrong, you want:

    x <- c(3 / 8, 1 / 2, 3 / 4, 1, 1 + 1 / 2, 2, 3, 6)  ## V1
    y <- c(2.4, 2.6, 2.8)  ## column names
    Z <- cbind(c(0.50, 0.59, 0.66, 0.71, 0.75, 0.78, 0.82, 0.87),  ## V2
               c(0.48, 0.57, 0.64, 0.69, 0.73, 0.76, 0.80, 0.85),  ## V3
               c(0.44, 0.53, 0.60, 0.65, 0.69, 0.72, 0.76, 0.81))  ## V4
    xp <- 3 / 4
    yp <- 2.70
    

    You already has a well defined matrix on a grid. For example, you can investigate your 3D data by:

    persp(x, y, Z)
    image(x, y, Z)
    contour(x, y, Z)
    

    I don't recommend pracma, as the interp2 function has a bug. I suggest interp.surface function from fields package for interpolation on a grid.

    library(fields)
    ## the list MUST has name `x`, `y`, `x`!
    ## i.e., unnamed list `list(x, y, Z)` does not work!
    interp.surface(list(x = x, y = y, z = Z), cbind(xp, yp))
    # [1] 0.62
    

    interp2 from pracma is inconsistent. The manual says Z matrix is length(x) by length(y), but the function really checks that Z must be length(y) by length(x).

    ## from manual
    
       Z: numeric ‘length(x)’-by-‘length(y)’ matrix.
    
    ## from source code of `interp2`
    
       lx <- length(x)
       ly <- length(y)
       if (ncol(Z) != lx || nrow(Z) != ly) 
           stop("Required: 'length(x) = ncol(Z)' and 'length(y) = nrow(Z)'.")
    

    So, in order to make interp2 works, you have to pass in the transpose of Z:

    interp2(x, y, t(Z), xp, yp)
    # [1] 0.62
    

    or reverse x and y (and xp, yp, too!!):

    interp2(y, x, Z, yp, xp)
    # [1] 0.62
    

    This is really inconsistent with the way that we work with image, contour and persp.