10

I have data in the form (x, y, z) where x and y are not on a regular grid. I wish to display a 2D colormap of these data, with intensity (say, grey scale) mapped to the z variable. An obvious solution is to interpolate (see below) on a regular grid,

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)

enter image description here

However, this loses the information of the initial mesh (position of the points with actual data), which could be very fine or very rough at certain locations. My preference would be for a delaunay tiling with triangles, which accurately represents the actual location and density of the original data points.

Ideally the solution would

  • compute the tesselation outside of the plotting function, so that the resulting polygons may be plotted with either ggplot2, lattice, or base graphics

  • be fast. In my real-life example (~1e5 points), the calculation of the tesselation via deldir can be really slow.

By "tesselation" I mean either Delaunay triangles or Voronoi diagrams, although my preference would be for the former. However it bring the additional complexity of interpolating the colour of each triangle based on the original data points.

baptiste
  • 75,767
  • 19
  • 198
  • 294

2 Answers2

4

Here's a solution based on dirichlet from the maptools package,

d <- data.frame(x=runif(1e3, 0, 30), y=runif(1e3, 0, 30))
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") 
plot(Z, col=grey(d$z/max(d$z)))

dirichlet

I'm still unsure of the way to extract the polygons from this SpatialPolygons class.

Also if there's an easy way to produce the "correct" colors for the associated delaunay tesselation I'd like to hear it.

baptiste
  • 75,767
  • 19
  • 198
  • 294
  • What do you mean by ""correct" colours"? Are you still trying to get the right grey shade attached to the right tile? If so, does the `@plotOrder` slot of Z help? By that I mean `plot(Z, col=grey(d$z/max(d$z))[Z@plotOrder])` – Gavin Simpson Apr 11 '11 at 07:47
  • Indeed, I tried this before I realised the real problem was that there are more delaunay tiles than z values. Try the above example with delaunay(), Y <- delaunay(X) Z <- as(Y, "SpatialPolygons") plot(Z, col=grey(d$z/max(d$z))[Z@plotOrder]) tiles get assigned the wrong color. I would need some kind of interpolation between adjacent tiles/points. – baptiste Apr 11 '11 at 22:05
2

Here is a lattice solution using deldir

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

pal1 <- grey(seq(0,1,leng=500))
library(latticeExtra)

 levelplot(z~x*y, data=d,
           panel = function(...) panel.voronoi(..., points=FALSE),
           interpolate=TRUE,
           col.regions = colorRampPalette(pal1)(1e3), cut=1e3)

enter image description here

baptiste
  • 75,767
  • 19
  • 198
  • 294