Interpolate irregular grid to regular grid

Question

I have an irregular grid, that I need to convert to a regular grid to take advantage of image's useRaster=TRUE option for graphics devices. I can do this on a small scale by converting the irregular grid to points, then interpolate the points using akima's interp. However this scales horribly with larger dimensions, so I'm looking for options.

First, here is the small scale (5x10) example, where only the x-dimension is irregular:

nx <- 5
ny <- 10
si <- list()  # irregular surface
si$x <- cumsum(runif(nx) * 10) + 100
si$y <- seq(20, 50, length.out=ny)
si$z <- matrix(rnorm(nx * ny), ncol=ny)
image(si)

And the bilinear interpolated result:

sr_x <- seq(min(si$x), max(si$x), length.out=nx * 5)
sr_y <- si$y  # this dimension is already regular
require(akima)  # interpolate from points repeated off irregular grid
sr <- interp(rep(si$x, length(si$y)), rep(si$y, each=length(si$x)), si$z,
             xo=sr_x, yo=sr_y)
image(sr, useRaster=TRUE)

However, if a larger dimension irregular grid is used (e.g. nx <- 50; ny <- 100), the procedure is really slow. Is there either a library or technique that would speed up the process?

---

Update and possibly a solution. The data describes time vs time (both in years), where the irregular dimension has timesteps between 0.5 days to 30 days, and the second time axis has equal-spaced 365 day spacings. Since the spacings are much smaller along the irregular axis, interpolating will not work. Thus a smoothing or aggregating method will yield better results.

A more realistic data scenario, showing the finer irregular dimension:

nx <- 200
ny <- 10
si <- list()  # irregular surface
si$x <- cumsum(runif(nx, 0.5, 30) / 365)
si$y <- 1:ny
si$z <- matrix(rnorm(nx * ny), ncol=ny)
image(si)

And some really crude aggregate means:

dx <- 1/12  # 1 month spacing along x-axis
sr <- list()  # regular surface
sr$x <- seq(min(si$x), max(si$y), dx)  # equal-width breaks
nsrx <- length(sr$x)
sr$y <- si$y  # this dimension is already regular
sr$z <- matrix(nrow=length(sr$x), ncol=length(sr$y))
# Classify irregular dimension
si_xc <- cut(si$x, sr$x, include.lowest=TRUE, labels=FALSE)
# Aggregate means from irregular to regular dimension
for(xi in seq_len(nsrx))
    sr$z[xi,] <- apply(si$z[si_xc == xi, , drop=FALSE], 2, mean)
image(sr, zlim=range(si$z), useRaster=TRUE)

This seems to do the trick, and scales on much larger datsets with 100s of years along each dimension. So I suppose my new question would be simply to tidy up the above code to perform aggregate means.

Answer 1

There are several packages w/ "kriging" tools, which is basically what you want. However, I dunno whether it'll be any faster than what akima::interp does.

I solved this using multicore techniques, so if you have a multicore processor, consider something similar to the following code snippet:

picbits <- clusterApply( myclus, 1:length(picsec) , function(j) { gc(); 
akima::interp(newx[picsec[[j]] ], newy[picsec[[j]] ], picture[picsec[[j]] ], 

xo=trunc(min(newx[picsec[[j]] ])):trunc(max(newx[picsec[[j]] ])), 

yo=trunc(min(newy[picsec[[j]] ])):trunc(max(newy[picsec[[j]] ])) )} )

That is extracted from a function I wrote to perform a rotational "swirl" on an image, so there's a lot of cruft there you won't need.