Thanks, it is cleaner indeed. However, it's still painfully slow on large data sets (typically a couple hundred thousand points). Looking at the source of various functions, it looks like the place to look at would be deldir::tile.list(). To my naive eyes it doesn't look particularly efficient, given that the "meat" of the calculation that I expected more costly was performed by deldir() in much less time,
N <- 1000 set.seed(123) d <- data.frame(x=rnorm(N, 15, 4), y=rnorm(N, 15, 4)) library(spatstat) library(maptools) W <- owin(c(0, 30), c(0, 30)) X <- as.ppp(d, W=W) system.time(dd <- deldir(d$x, d$y) ) user system elapsed 0.079 0.009 0.088 system.time(Y <- dirichlet(X)) user system elapsed 5.494 0.074 5.756 tile.list() would probably benefit greatly from being ported to C++, and currently returns a list where I really want a data.frame to avoid subsequent looping. Further coercion steps that I attempted earlier are unnecessary as I wouldn't use the tess or SpatialPolygons classes. Suggestions and comments most welcome. baptiste On 13 April 2011 14:36, Hadley Wickham <[email protected]> wrote: > Did you try fortify on the SpatialPolygons? > Hadley > > On Tue, Apr 12, 2011 at 4:24 PM, baptiste auguie > <[email protected]> wrote: >> 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 >> >> _______________________________________________ >> R-sig-Geo mailing list >> [email protected] >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > > > -- > Assistant Professor / Dobelman Family Junior Chair > Department of Statistics / Rice University > http://had.co.nz/ > _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
