"interp" for discrete points to get heatmap/contour R

Question

I have a sitation in which I generate data from simulation and then would like to plot (heat map/contour/3d plot etc); however, for these, data needs to be interpolated using functions like interp. Here is the sample dataset.

Here is the piece of code I tried...

library(akima)
library(GA) # for persp3D; there exists another package for same function "fields"

data <- read.table(commandArgs()[3], header=T,sep="\t")
data <- na.omit(data)

qmax = max(data$q)
kmax = max(data$k)

x <- data$k_bike/kmax
y <- data$k_car/kmax
z <- data$q/qmax

matrix = interp(x,y,z)

persp3D(matrix ,nlevels=30, asp=1, xlim=c(0,1), ylim=c(0,1), color.palette=colorRampPalette(c("green3","yellow", "red"),space = "rgb") )

so the result is --

Now, due to interpolation, there are many points, which have red/orange color instead of green or so. For e.g, if I use levelplot of lattice

levelplot(z~x*y, xlim=c(0,1), ylim=c(0,1), col.regions=colorRampPalette(c("green3","yellow", "red"),space = "rgb") )

The outcome is --

Now, it is clearly visible that there are very few data points having zero (or almost zero) zvalue. Now, the problem is, with levelplot, I get artefacts (white color for missing data points) and I would like to have a better interpolation. Is there any other function to perform this?

I also tried contour plots as follows:

  scale <- (qmax+10) / qmax * c(0.000, 0.01, 0.05, 0.10, 0.25, 0.5, 0.75, 1.0)
filled.contour(matrix, nlevels=30, asp=1, xlim=c(0,1), ylim=c(0,1), levels=scale,color.palette=colorRampPalette(c("green3","yellow", "red"),space = "rgb") )

and result is again (kind of wrong color indication).

In short -- I would like to have contour plot or 3d plot but with a clear (or correct) color indication of zero (about zero) zvalue data points similar to level plot.

Answer 1

I approached your question with deldir and rgl packages (they allow plotting of surfaces defined by irregular collections of points).

library(deldir); library(rgl)

# Below two lines require time depending on the　machine power, be careful
dxyz <- deldir(x, y, z = z)      # do Delaunay triangulation
mxyz <- as.mesh3d(dxyz)          # convert it to triangle.mesh3d.obj

bgyr <- colorRampPalette(c("blue", "green", "yellow", "red"))    # colour func

# use z values for colouring
plot3d(mxyz, col=bgyr(256)[cut(mxyz$vb[3,], 256)][mxyz$it], type="shade")
light3d()          # if you want vivit colours

# another approach
# you can solve it by just increasing interp()'s arguments, nx and ny.

library(akima); library(lattice); library(dplyr)

df <- interp(x,y,z, nx=150, ny=150) %>% interp2xyz() %>% data.frame()
levelplot(z ~ x * y, df, xlim=c(0,1), ylim=c(0,1), 
          col.regions = colorRampPalette(c("green3", "yellow",  "red"), space = "rgb"))