3D Surface Interpolation

Question

I have a three column data frame with latitude, longitude, and underground measurements as the columns. I am trying to figure out how to interpolate data points between the points I have (which are irregularly space) and then create a smooth surface plot of the entire area. I have tried to use the 'surface3d' function in the 'rgl' package but my result looks like a single giant spike. I have been able to plot the data with 'plot3d' but I need to take it a step further and fill in the blank spaces with interpolation. Any ideas or suggestions? I'm also open to using other packages, the rgl just seemed like the best fit at the time.

EDIT: here's an excerpt from my data (measurements of aquifer depth) :

 lat_dd_NAD83 long_dd_NAD83 lev_va_ft
1     37.01030     -101.5006    288.49
2     37.03977     -101.6633    191.68
3     37.05201     -100.4994    159.34
4     37.06567     -101.3292    174.07
5     37.06947     -101.4561    285.08
6     37.10098     -102.0134    128.94

You can use the `loess.surf` function in the `asbio` package. If you could provide your data we can help you better — adaien, Apr 04 '16 at 22:09
kriging (with gstat package) would seem an obvious technique to try — Richard Telford, Apr 04 '16 at 22:39
Various interpretations of "3D surface" are possible. Are you trying to interpolate between the measured `lev` values on the basis of the `lat` and `long`? (I would call that a 2d manifold embedded in 3 dimensions. And if so then the suggestion of loess (albeit in the stats package) with 2 variables on the RHS seems to be on point and has already an illustrated answer: http://stackoverflow.com/questions/15019725/r-3d-surface-plot — IRTFM, Apr 05 '16 at 05:51

Dmitriy · Answer 1 · 2016-11-25T16:14:56.327

Just to add small but (maybe) important note about interpolation.

Using very nice package "akima" you can easily interpolate your data:

library(akima)
library(rgl)
# library(deldir)

# Create some fake data
x <- rnorm(100)
y <- rnorm(100)
z <- x^2 + y^2

# # Triangulate it in x and y
# del <- deldir(x, y, z = z)
# triangs <- do.call(rbind, triang.list(del))
# 
# # Plot the resulting surface
# plot3d(x, y, z, type = "n")
# triangles3d(triangs[, c("x", "y", "z")], col = "gray")

n_interpolation <- 200

spline_interpolated <- interp(x, y, z,
                              xo=seq(min(x), max(x), length = n_interpolation),
                              yo=seq(min(y), max(y), length = n_interpolation),
                              linear = FALSE, extrap = TRUE)

x.si <- spline_interpolated$x
y.si <- spline_interpolated$y
z.si <- spline_interpolated$z

persp3d(x.si, y.si, z.si, col = "gray")

Spline - interpolated picture (200 steps)

With this package you can easily change amount of steps of interpolation, etc. You will need at least 10 (the more the better) points to get a reasonable spline interpolation with this package. Linear version works well regardless amount of points.

P.S. Thanks for user 2554330 - didn't knew about deldir, really useful thing in some cases.

user2554330 · Accepted Answer · 2016-04-13T12:06:44.033

You could use the deldir package to get a Delaunay triangulation of your points, then convert it to the form of data required by triangles3d for plotting. I don't know how effective this would be on a really large dataset, but it seems to work on 100 points:

library(deldir)
library(rgl)
# Create some fake data
x <- rnorm(100)
y <- rnorm(100)
z <- x^2 + y^2

# Triangulate it in x and y
del <- deldir(x, y, z = z)
triangs <- do.call(rbind, triang.list(del))

# Plot the resulting surface
plot3d(x, y, z, type = "n")
triangles3d(triangs[, c("x", "y", "z")], col = "gray")

EDITED to add:

The version of rgl on R-forge now has a function to make this easy. You can now produce a plot similar to the one above using

library(deldir)
library(rgl)
plot3d(deldir(x, y, z = z))

There is also a function to construct mesh3d objects from the deldir() output.

3D Surface Interpolation

2 Answers2