Tenho dados no formato (x, y, z) em que x e y não estão em uma grade regular. Desejo exibir um mapa de cores 2D desses dados, com intensidade (digamos, escala de cinza) mapeada para a variável z. Uma solução óbvia é interpolar (veja abaixo) em uma grade 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)

No entanto, isso perde as informações da malha inicial (posição dos pontos com dados reais), que podem ser muito finas ou muito irregulares em determinados locais. Minha preferência seria por um mosaico de delaunay com triângulos, que representa com precisão a localização real e a densidade dos pontos de dados originai

Ideally a solução seria

calcule o mosaico fora da função de plotagem, para que os polígonos resultantes possam ser plotados comggplot2, lattice ou gráficos de base

seja rápido. No meu exemplo da vida real (~ 1e5 pontos), o cálculo do mosaico viadeldir pode ser bem lento.

Por "mosaico", quero dizer triângulos de Delaunay ou diagramas de Voronoi, embora minha preferência seja pelo primeiro. No entanto, traz a complexidade adicional de interpolar a cor de cada triângulo com base nos pontos de dados originai