R/streamplot.R

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
              }
              ## 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)
              }
            ## 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)
            }

            ## 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
          }
          )


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,
            ...){
              if (missing(layers)) layers=seq_len(nlayers(object))

              if (isField == 'dXY'){
                  isField <- TRUE
                  u <- subset(object, 1)
                  v <- subset(object, 2)
                  slope <- sqrt(u^2 + v^2)
                  aspect <- atan2(v, u)
                  aspect <- (pi/2 - aspect) %% (2 * pi)
                  object <- stack(slope, aspect)
              }
              
              if (isTRUE(isField))
                  callNextMethod(object, layers, droplet,
                                 streamlet,
                                 par.settings,
                                 isField=TRUE,
                                 reverse,
                                 parallel,
                                 mc.cores, cl, ...)
            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,
                               isField=FALSE, reverse=reverse,
                               parallel=parallel, mc.cores=mc.cores,
                               cl, ...)
              p
          }
          })

Try the rasterVis package in your browser

Any scripts or data that you put into this service are public.

rasterVis documentation built on May 29, 2017, 1:05 p.m.