Skip to content

Commit

Permalink
Allow for the colorkey in streamplot as proposed in #38
Browse files Browse the repository at this point in the history
  • Loading branch information
@oscarperpinan
oscarperpinan committed Dec 23, 2016
1 parent d0aa0c0 commit ea1c49e
Show file tree
Hide file tree
Showing 3 changed files with 186 additions and 175 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: rasterVis
Type: Package
Title: Visualization Methods for Raster Data
Version: 0.41
Date: 2016-12-10
Version: 0.42
Date: 2016-12-23
Encoding: UTF-8
Authors@R: c(person("Oscar", "Perpinan Lamigueiro", email="oscar.perpinan@gmail.com", role=c('cre', 'aut')), person("Robert", "Hijmans", email= "r.hijmans@gmail.com", role='aut'))
Description: Methods for enhanced visualization and interaction with raster data. It implements visualization methods for quantitative data and categorical data, both for univariate and multivariate rasters. It also provides methods to display spatiotemporal rasters, and vector fields. See the website for examples.
Expand Down
340 changes: 172 additions & 168 deletions R/streamplot.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,172 +3,174 @@ setGeneric('streamplot', function(object, ...){standardGeneric('streamplot')})
setMethod('streamplot',
signature(object='Raster'),
definition = function(object, layers,
droplet = list(), streamlet = list(),
par.settings=streamTheme(),
isField = FALSE, reverse=FALSE, ##unit = 'radians',
parallel=TRUE, mc.cores=detectCores(), cl=NULL,
...){
stopifnot(is.list(droplet))
stopifnot(is.list(streamlet))
## uniVector extracts the unitary vector at point
## xy. It needs xrange, yrange and aspectVals as
## "global" variables
uniVector <- function(xy){
idx <- findInterval(xy[1], xrange, rightmost.closed=TRUE, all.inside=TRUE)
idy <- findInterval(xy[2], rev(yrange), rightmost.closed=TRUE, all.inside=TRUE)
idy <- nrow(aspectVals) - idy + 1
ang <- aspectVals[idy, idx]
c(sin(ang), cos(ang))
}
## rk4 calculates the integral curve with Runge-Kutta
## order 4 from point p with step h
rk4 <- function(p, h){
xy <- p
v <- uniVector(p)
k1 <- h*v
v2 <- uniVector(xy + 1/2*k1)
k2 <- h*v2
v3 <- uniVector(xy + 1/2*k2)
k3 <- h*v3
v4 <- uniVector(xy + k3)
k4 <- h*v4
xy + 1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4
}
## streamLine creates a streamlet starting from point
## p, with L steps of length h both in forward and
## backward directions. Starting point is included in
## the result, therefore total length is 2*(L*h) + 1).
streamLine <- function(p, h=0.5, L=10, pal, cex){
## Color class
colClass <- p[3]
## Coordinates
p <- matrix(p[1:2], ncol=2)
## Forward direction
forw <- matrix(nrow=L, ncol=2)
xy <- p
for (i in 1:L){
xy <- rk4(xy, h=h)
forw[i,] <- xy
droplet = list(), streamlet = list(),
par.settings=streamTheme(),
colorkey = FALSE,
isField = FALSE, reverse=FALSE, ##unit = 'radians',
parallel=TRUE, mc.cores=detectCores(), cl=NULL,
...){
stopifnot(is.list(droplet))
stopifnot(is.list(streamlet))
## uniVector extracts the unitary vector at point
## xy. It needs xrange, yrange and aspectVals as
## "global" variables
uniVector <- function(xy){
idx <- findInterval(xy[1], xrange, rightmost.closed=TRUE, all.inside=TRUE)
idy <- findInterval(xy[2], rev(yrange), rightmost.closed=TRUE, all.inside=TRUE)
idy <- nrow(aspectVals) - idy + 1
ang <- aspectVals[idy, idx]
c(sin(ang), cos(ang))
}
## Backward direction
back <- matrix(nrow=L, ncol=2)
xy <- p
for (i in 1:L){
xy <- rk4(xy, h=-h)
back[i,] <- xy
## rk4 calculates the integral curve with Runge-Kutta
## order 4 from point p with step h
rk4 <- function(p, h){
xy <- p
v <- uniVector(p)
k1 <- h*v
v2 <- uniVector(xy + 1/2*k1)
k2 <- h*v2
v3 <- uniVector(xy + 1/2*k2)
k3 <- h*v3
v4 <- uniVector(xy + k3)
k4 <- h*v4
xy + 1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4
}
## Streamlet: backward -> starting point -> forward
pts <- rbind(back[L:1,], p, forw)
## Only non-NA points are useful
whichNA <- apply(pts, 1, function(x)any(is.na(x)))
pts <- pts[!whichNA,]


## Perhaps only the starting point remains...
if (!is.matrix(pts)) stream <- list(x=pts[1], y=pts[2])
else stream <- list(x=pts[,1], y=pts[,2])
## Color and size
ncolors <- length(stream$x)
color <- pal[colClass]
stream$colors <- rev(colorRampPalette(c(color, 'gray50'))(ncolors))
stream$cexs <- seq(.3, cex, length=ncolors)
stream
}

if (isField) {
field <- object
## streamLine creates a streamlet starting from point
## p, with L steps of length h both in forward and
## backward directions. Starting point is included in
## the result, therefore total length is 2*(L*h) + 1).
streamLine <- function(p, h=0.5, L=10, pal, cex){
## Color class
colClass <- p[3]
## Coordinates
p <- matrix(p[1:2], ncol=2)
## Forward direction
forw <- matrix(nrow=L, ncol=2)
xy <- p
for (i in 1:L){
xy <- rk4(xy, h=h)
forw[i,] <- xy
}
## Backward direction
back <- matrix(nrow=L, ncol=2)
xy <- p
for (i in 1:L){
xy <- rk4(xy, h=-h)
back[i,] <- xy
}
## Streamlet: backward -> starting point -> forward
pts <- rbind(back[L:1,], p, forw)
## Only non-NA points are useful
whichNA <- apply(pts, 1, function(x)any(is.na(x)))
pts <- pts[!whichNA,]


## Perhaps only the starting point remains...
if (!is.matrix(pts)) stream <- list(x=pts[1], y=pts[2])
else stream <- list(x=pts[,1], y=pts[,2])
## Color and size
ncolors <- length(stream$x)
color <- pal[colClass]
stream$colors <- rev(colorRampPalette(c(color, 'gray50'))(ncolors))
stream$cexs <- seq(.3, cex, length=ncolors)
stream
}

if (isField) {
field <- object
} else {
field <- terrain(object, opt = c("slope", "aspect"))##, unit=unit)
field <- terrain(object, opt = c("slope", "aspect"))##, unit=unit)
}
## Values to be used by uniVector
xrange <- xFromCol(field, 1:ncol(field))
yrange <- yFromRow(field, 1:nrow(field))
aspectVals <- as.matrix(subset(field, 2))
## Should streamlets go from sinks to sources?
if (reverse) aspectVals <- aspectVals + pi

## Droplets and streamlets configuration
default.droplet <- list(cropExtent = .97, pc = .5)
droplet <- modifyList(default.droplet, droplet)

default.streamlet <- list(L=10, h = mean(res(object)))
streamlet <- modifyList(default.streamlet, streamlet)

## Crop original object to avoid droplets at boundaries
cropField <- crop(field, extent(field)*droplet$cropExtent)
nDroplets <- ncell(cropField) * droplet$pc/100
## Build a regular grid ...
texture <- sampleRegular(cropField, nDroplets, sp=TRUE)
## but only with coordinates where slope (texture[[1]]
## is not NA
crds <- coordinates(texture[!is.na(texture[[1]]),])
## Add jitter to regular grid
crdsNoise <- cbind(jitter(crds[,1]), jitter(crds[,2]))
## "Grid" of droplets
texture <- SpatialPointsDataFrame(crdsNoise,
data.frame(extract(field, crdsNoise)))
pts <- data.frame(t(coordinates(texture)))
npts <- length(texture)

## Color classes
pars <- if (is.function(par.settings)) {
par.settings()$superpose.symbol
} else {
par.settings$superpose.symbol
}
slopeVals <- texture[[1]]
nClasses <- length(pars$col)
colClasses <- pretty(slopeVals, n=nClasses)
indCol <- findInterval(slopeVals, colClasses,
rightmost.closed=TRUE, all.inside=TRUE)

## Stream Lines Calculation
pts <- rbind(pts, indCol)
h <- streamlet$h
L <- streamlet$L

## Forking does not work with Windows
if (.Platform$OS.type == "windows" & is.null(cl)) parallel <- FALSE

if (isTRUE(parallel)) {
if (!is.null(cl)) { ## parallel with a cluster
streamList <- parLapply(cl, pts, streamLine, h=h, L=L,
pal=pars$col, cex=pars$cex)
} else { ## parallel with forking
streamList <- mclapply(pts, streamLine, h=h, L=L,
pal=pars$col, cex=pars$cex,
mc.cores=mc.cores)
## Values to be used by uniVector
xrange <- xFromCol(field, 1:ncol(field))
yrange <- yFromRow(field, 1:nrow(field))
aspectVals <- as.matrix(subset(field, 2))
## Should streamlets go from sinks to sources?
if (reverse) aspectVals <- aspectVals + pi

## Droplets and streamlets configuration
default.droplet <- list(cropExtent = .97, pc = .5)
droplet <- modifyList(default.droplet, droplet)

default.streamlet <- list(L=10, h = mean(res(object)))
streamlet <- modifyList(default.streamlet, streamlet)

## Crop original object to avoid droplets at boundaries
cropField <- crop(field, extent(field)*droplet$cropExtent)
nDroplets <- ncell(cropField) * droplet$pc/100
## Build a regular grid ...
texture <- sampleRegular(cropField, nDroplets, sp=TRUE)
## but only with coordinates where slope (texture[[1]]
## is not NA
crds <- coordinates(texture[!is.na(texture[[1]]),])
## Add jitter to regular grid
crdsNoise <- cbind(jitter(crds[,1]), jitter(crds[,2]))
## "Grid" of droplets
texture <- SpatialPointsDataFrame(crdsNoise,
data.frame(extract(field, crdsNoise)))
pts <- data.frame(t(coordinates(texture)))
npts <- length(texture)

## Color classes
pars <- if (is.function(par.settings)) {
par.settings()$superpose.symbol
} else {
par.settings$superpose.symbol
}
slopeVals <- texture[[1]]
nClasses <- length(pars$col)
colClasses <- pretty(slopeVals, n=nClasses)
indCol <- findInterval(slopeVals, colClasses,
rightmost.closed=TRUE, all.inside=TRUE)

## Stream Lines Calculation
pts <- rbind(pts, indCol)
h <- streamlet$h
L <- streamlet$L

## Forking does not work with Windows
if (.Platform$OS.type == "windows" & is.null(cl)) parallel <- FALSE

if (isTRUE(parallel)) {
if (!is.null(cl)) { ## parallel with a cluster
streamList <- parLapply(cl, pts, streamLine, h=h, L=L,
pal=pars$col, cex=pars$cex)
} else { ## parallel with forking
streamList <- mclapply(pts, streamLine, h=h, L=L,
pal=pars$col, cex=pars$cex,
mc.cores=mc.cores)
}
} else { ## without parallel
streamList <- lapply(pts, streamLine, h=h, L=L,
pal=pars$col, cex=pars$cex)
}
} else { ## without parallel
streamList <- lapply(pts, streamLine, h=h, L=L,
pal=pars$col, cex=pars$cex)
}

## Points with low slope values are displayed earlier
idOrdered <- order(slopeVals)
streamList <- streamList[order(slopeVals)]

p <- levelplot(object, layers=1,
margin=FALSE, colorkey=FALSE,
par.settings=par.settings, ...) +
layer(lapply(streamList, function(streamlet){
panel.points(streamlet$x, streamlet$y,
col=streamlet$colors, cex=streamlet$cexs)
}), data= list(streamList=streamList))

p

## Points with low slope values are displayed earlier
idOrdered <- order(slopeVals)
streamList <- streamList[order(slopeVals)]

p <- levelplot(object, layers = 1,
margin = FALSE, colorkey = colorkey,
par.settings=par.settings, ...) +
layer(lapply(streamList, function(streamlet){
panel.points(streamlet$x, streamlet$y,
col=streamlet$colors, cex=streamlet$cexs)
}), data= list(streamList=streamList))

p
}
)


setMethod('streamplot',
signature(object='RasterStack'),
definition = function(object, layers,
droplet = list(), streamlet = list(),
par.settings=streamTheme(),
isField = FALSE, reverse=FALSE, ##unit = 'radians',
parallel=TRUE, mc.cores=detectCores(), cl=NULL,
...){
droplet = list(), streamlet = list(),
par.settings=streamTheme(),
colorkey = FALSE,
isField = FALSE, reverse=FALSE, ##unit = 'radians',
parallel=TRUE, mc.cores=detectCores(), cl=NULL,
...){
if (missing(layers)) layers=seq_len(nlayers(object))

if (isField == 'dXY'){
Expand All @@ -185,22 +187,24 @@ setMethod('streamplot',
callNextMethod(object, layers, droplet,
streamlet,
par.settings,
isField=TRUE,
colorkey,
isField = TRUE,
reverse,
parallel,
mc.cores, cl, ...)
else {
if (!missing(layers)) {
object <- subset(object, subset=layers)
else {
if (!missing(layers)) {
object <- subset(object, subset=layers)
}
objectList <- unstack(object)
names(objectList) <- names(object)
p <- xyplot.list(objectList, FUN = streamplot,
droplet = droplet, streamlet = streamlet,
par.settings = par.settings,
colorkey = colorkey,
isField = FALSE, reverse = reverse,
parallel = parallel, mc.cores = mc.cores,
cl, ...)
p
}
objectList <- unstack(object)
names(objectList) <- names(object)
p <- xyplot.list(objectList, FUN=streamplot,
droplet=droplet, streamlet=streamlet,
par.settings=par.settings,
isField=FALSE, reverse=reverse,
parallel=parallel, mc.cores=mc.cores,
cl, ...)
p
}
})
Loading

0 comments on commit ea1c49e

Please sign in to comment.