# Interpolating non-gridded data

**Saturn Elephant**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I still wrapped a part of the C++ library **CGAL** in a R
package, namely
**interpolation**.

The purpose of this package is to perform interpolation of bivariate functions. As compared to existing packages, it can do more: it can interpolate vector-valued functions (with dimension two or three), and it does not require that the given data are gridded. I will illustrate this second point here.

First, let’s plot a surface \(z = f(x, y)\).

# make a square grid x <- y <- seq(-3, 3, by = 0.1) Grid <- expand.grid(X = x, Y = y) # make surface z <- with(Grid, 30 * dcauchy(X) * dcauchy(Y)) # plot library(deldir) delxyz <- deldir(Grid[["X"]], Grid[["Y"]], z = z) library(rgl) open3d(windowRect = 50 + c(0, 0, 512, 512)) persp3d(delxyz, color = "red")

Now we will make a hole in this surface, and then we will interpolate it.

# make a hole in the grid: the square [-0.5, 1]x[-0.5, 1] toremove <- with(Grid, X > -0.5 & X < 1 & Y > -0.5 & Y < 1) GridWithHole <- Grid[!toremove, ] # plot this grid plot(GridWithHole[["X"]], GridWithHole[["Y"]], pch = 19, asp = 1)

Now, to plot the surface with the hole, I will use a constrained Delaunay triangulation. I didn’t find a more straightforward way.

# constraint edges: squares [-0.5, 1]x[-0.5, 1] and [-3, 3]x[-3, 3] Constraints <- data.frame( X = c(-0.5, 1, 1, -0.5, -3, 3, 3, -3), Y = c(-0.5, -0.5, 1, 1, -3, -3, 3, 3) ) # add constraint edges to the grid with the hole GridWithHole <- rbind(Constraints, GridWithHole) # remove duplicated points GridWithHole <- GridWithHole[!duplicated(GridWithHole), ] # vertices points <- cbind(GridWithHole[["X"]], GridWithHole[["Y"]]) # constraint edges must be given by vertex indices edges <- rbind(c(1L, 2L), c(2L, 3L), c(3L, 4L), c(4L, 1L)) # the hole edges <- rbind(edges, edges + 4L) # the outer square # perform constrained Delaunay triangulation library(delaunay) del <- delaunay(points, constraints = edges)

Note that the **delaunay** package is also a wrapper of
**CGAL**.

This Delaunay triangulation provides triangular faces that we can use to
create a 3D **rgl** mesh.

# create surface plot z <- with(GridWithHole, 30 * dcauchy(X) * dcauchy(Y)) vertices <- cbind(points, z) rmesh <- tmesh3d( vertices = t(vertices), indices = t(del[["faces"]]), homogeneous = FALSE ) rmesh <- addNormals(rmesh) # we plot the front side in red and the back side in gray open3d(windowRect = 50 + c(0, 0, 512, 512)) shade3d(rmesh, color = "red", back = "cull") shade3d(rmesh, color = "gray", front = "cull") persp3d(delxyz, add = TRUE, alpha = 0.2)

Good. Now let’s interpolate.

# make the interpolating function library(interpolation) fun <- interpfun( GridWithHole[["X"]], GridWithHole[["Y"]], z, method = "linear" ) # make a grid of the hole xnew <- ynew <- seq(-0.5, 1, by = 0.04) griddedHole <- expand.grid(X = xnew, Y = ynew) # interpolation znew <- fun(griddedHole[["X"]], griddedHole[["Y"]]) # new 3D data newPoints <- cbind(griddedHole[["X"]], griddedHole[["Y"]], znew) # plot open3d(windowRect = 50 + c(0, 0, 512, 512)) shade3d(rmesh, color = "red", back = "cull") shade3d(rmesh, color = "gray", front = "cull") points3d(newPoints)

Not very nice, you think? Right, but I used the linear method of
interpolation here. The **interpolation** package also
provides the *Sibson* method, this one is not linear. One just
has to repeat the above code but starting with:

fun <- interpfun( GridWithHole[["X"]], GridWithHole[["Y"]], z, method = "sibson" )

And we obtain:

This is not exactly the true curve, but nevertheless this is impressive.

**leave a comment**for the author, please follow the link and comment on their blog:

**Saturn Elephant**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.