[R-sig-Geo] colour mapping tesselated irregular xy data

[baptiste auguie]

Dear list,

I have x,y,z data originating from an electromagnetic simulation on an
irregular mesh. I wish to visualise a 2D map of these (N~1e6) data
points with colour-coded z intensity.

On way is to perform an interpolation onto a regular grid. This works
(using akima's interp for example), but is not what I'm after. The xy
mesh is very fine at places, and rough at others, and interpolation
loses this spatial information. Thus, I wish to display coloured tiles
of irregular size. Voronoi and delaunay tesselations seem quite suited
to the task; however I am stuck with two technical difficulties.

1- while voronoi/dirichlet gives me N polygons, delaunay obviously
returns more triangular tiles than there are points initially. How
could I perform the interpolation for the assignment of colours?

2- spatstat and maptools seem to work with special|spatial| classes,
from which I have been unable to efficiently extract the polygon
outlines in a data.frame. This is necessary because I wish to use
lattice or ggplot2 for plotting, rather than the default plot method.

Below's a minimal illustration,

N <- 100
set.seed(123)
d <- data.frame(x=rnorm(N, 15, 4), y=rnorm(N, 15, 4))
## a circular intensity pattern
d$z = (d$x - 15)^2 + (d$y - 15)^2

library(spatstat)
library(maptools)

W <- ripras(df, shape="rectangle")
W <- owin(c(0, 30), c(0, 30))
X <- as.ppp(d, W=W)
Y <- dirichlet(X)
Z <- as(Y, "SpatialPolygons")
## default plot
## plot(Z, col=grey(d$z/max(d$z)))

## awfully inefficient extraction of the polygons...

extract.tiles <- function(x){
 data.frame(x$bdry[[1]][1:2])
}

all <- lapply(Y$tiles, extract.tiles )

library(ggplot2)
## convert list to long format data.frame with ids
names(all) <- seq_along(all)
m <- melt(all, id=1:2)
m$fill <- d$z[as.integer(m$L1)]

ggplot(m)+
 geom_polygon(aes(x,y,fill=fill, group=fill))+
 scale_fill_gradient(low="black", high="white")

library(latticeExtra)
## result similar (except for the edges) to latticeExtra's
panel.voronoi, as intended
levelplot(z~x*y, data=d,
         panel = function(...) panel.voronoi(..., points=FALSE),
         col.regions = colorRampPalette(grey(0:1))(1e3), cut=1e3)


## now trying the same with delaunay tiles
Y1 <- delaunay(X)
Z1 <- as(Y1, "SpatialPolygons")
## no link between original z-values and tiles
plot(Z1, col=grey(d$z/max(d$z)))


Best regards,

baptiste