[R] Filled vector contours

Jan Tosovsky j.tosovsky at email.cz
Mon Aug 25 22:17:44 CEST 2014


On 2014-08-25 David Winsemius wrote:
> On 2014-08-24 Jan Tosovsky wrote:
> > On 2014-08-24 Prof Brian Ripley wrote:
> > > On 24/08/2014 08:51, Jan Tosovsky wrote:
> > > >
> > > > I am trying to create vector output (SVG) of filled contours.
> > > >
> > > > I've found two approaches so far:
> > > >
> > > > (1)
> > > > library('lattice')
> > > > svg("D:/test.svg")
> > > > filled.contour(volcano)
> > > > #levelplot(volcano, panel=panel.levelplot.raster)
> > > > dev.off()
> > > >
> > > > (2)
> > > > library("raster")
> > > > svg("D:/test.svg")
> > > > rr <- raster(t(volcano))
> > > > rc <- cut(rr, breaks= 10)
> > > > pols <- rasterToPolygons(rc, dissolve=T)
> > > > spplot(pols)
> > > > dev.off()
> > > >
> > > > But I'd like to get smooth closed polygons not broken into small
> > > > cells.
> > >
> > > But the region between two contours is not in general a closed
> > > polygon.
> > 
> > I appologize for not being precise here.
> > 
> > My goal is to get something like this:
> > http://drifted.in/other/contours/composition.svg
> > 
> > which is a composition of level plots, e.g.
> > http://drifted.in/other/contours/level_plot.svg
> > 
> > By 'level plot' I mean the cut of the data at certain level. It is
> > always a closed polygon (or multiple polygons) consisting of contour 
> > lines connected, if required, by border lines.
> > 
> > For me it is preferred way over 'isoband', which doesn't include areas
> > hidden by upper levels:
> > http://drifted.in/other/contours/isoband.svg
> > 
> 
> Most people with any R experience would have thought that "levelplot"
> referred to something like:
> (Example 6.9 from Lattice, by Deepayan Sarkar to which I have only
> added the contour line parameter which it appears is what you were 
> seeking.)
>
> library("locfit") ; library(lattice)
> env <- environmental
> env$ozone <- env$ozone^(1/3)
> env$Radiation <- equal.count(env$radiation, 4)
> fm1.env <- lm(ozone ~ radiation * temperature * wind, env)
> fm2.env <- loess(ozone ~ wind * temperature * radiation, env, span =
> 0.75, degree = 1)
> fm3.env <- loess(ozone ~ wind * temperature * radiation, env,
> parametric
> = c("radiation", "wind"), span = 0.75, degree = 2)
> > fm4.env <- locfit(ozone ~ wind * temperature * radiation, env)
> w.mesh <- with(env, do.breaks(range(wind), 50))
> t.mesh <- with(env, do.breaks(range(temperature), 50))
> r.mesh <- with(env, do.breaks(range(radiation), 3))
> grid <- expand.grid(wind = w.mesh, temperature = t.mesh, radiation =
> r.mesh)
> grid[["fit.linear"]] <- predict(fm1.env, newdata = grid)
> grid[["fit.loess.1"]] <- as.vector(predict(fm2.env, newdata = grid))
> grid[["fit.loess.2"]] <- as.vector(predict(fm3.env, newdata = grid))
> grid[["fit.locfit"]] <- predict(fm4.env, newdata = grid)
> 
> png()
> print(levelplot(fit.linear + fit.loess.1 + fit.loess.2 + fit.locfit ~
> wind * temperature | radiation, data = grid, contour=TRUE))
> dev.off()

Thanks, it looks similar, but in reality countours are just polylines not
connected into polygons and the fill is made from individual cells. But good
to know there is another way. I've realized it is a variant of commented
line in (1). So here is a simplified version:

(3)
library('lattice')
svg("D:/test.svg")
levelplot(volcano, contour=TRUE)
dev.off()

In the meantime I've found fourth one (based on Grid2Polygons library
example):

library(sp)
library(Grid2Polygons)
data(meuse.grid)
coordinates(meuse.grid) <- ~ x + y
gridded(meuse.grid) <- TRUE
meuse.grid <- as(meuse.grid, "SpatialGridDataFrame")
meuse.plys.lev <- Grid2Polygons(meuse.grid, "dist", level = TRUE)

svg("D:/meuse.svg")
plot(meuse.plys.lev, col = heat.colors(length(meuse.plys.lev)))
dev.off()

It is similar to (2). It actually doesn't produce smooth countours, just
joins cells into areas of the same fill/level. And these areas are
'isobands', not 'level plots' (or rather isolevels?).

I am leaving this open. For the future reference it would be nice to modify
code (4) for volcano data. I am lost here ;-)

Thanks, Jan



More information about the R-help mailing list