plotar y colorear datos en cuadrícula irregular

Tengo datos en la forma (x, y, z) donde x e y no están en una cuadrícula regular. Deseo mostrar un mapa de color 2D de estos datos, con intensidad (por ejemplo, escala de grises) asignada a la variable z. Una solución obvia es interpolar (ver más abajo) en una cuadrícula regular,

d <- data.frame(x=runif(1e3, 0, 30), y=runif(1e3, 0, 30))
d$z = (d$x - 15)^2 + (d$y - 15)^2


library(akima)
d2 <- with(d, interp(x, y, z, xo=seq(0, 30, length = 30),
                     yo=seq(0, 30, length = 50), duplicate="mean"))

pal1 <- grey(seq(0,1,leng=500))
with(d2, image(sort(x), sort(y), z, useRaster=TRUE, col = pal1))
points(d$x, d$y, col="white", bg=grey(d$z/max(d$z)), pch=21, cex=1,lwd=0.1)

Sin embargo, esto pierde la información de la malla inicial (posición de los puntos con datos reales), que podría ser muy fina o muy rugosa en ciertos lugares. Mi preferencia sería un mosaico delaunay con triángulos, que representa con precisión la ubicación real y la densidad de los puntos de datos originales.

Idealmente la solución sería

computar la teselación fuera de la función de trazado, de modo que los polígonos resultantes se puedan trazar conggplot2, lattice, o gráficos base

e rápido. En mi ejemplo de la vida real (~ 1e5 puntos), el cálculo de la teselación a través dedeldir puede ser muy lento.

Por "teselación" me refiero a triángulos de Delaunay o diagramas de Voronoi, aunque mi preferencia sería la primera. Sin embargo, aporta la complejidad adicional de interpolar el color de cada triángulo en función de los puntos de datos originales.

Respuestas a la pregunta(2)

Su respuesta a la pregunta