[R-sig-Geo] raster (graphics) support in R 2.11.0

[Michael Sumner]

Hello, I've been experimenting with the new raster graphics support in
2.11.0 (unstable)  It's nice as it can provide graphical
interpolation, and arbitrary rotation (angle is not obviously useful
in this context, but I've included it for experimenting).

This is a version of sp:::image.SpatialGridDataFrame that sets up the
data for the raster graphics function. I'm not sure how to handle
missing values, here they are just set to white. Probably it's best to
live with the warnings from raster, but it might be possible to
determine if the device can deal with it. "Per-pixel transparency" is
not supported on the windows() or png() device (on Windows) but it
seems to work for pdf().

The raster() function is very raw so it takes a bit of setting up to
use. It expects an existing plot that is added to, and it can take a
matrix of (hex) colours or a matrix or 3D array of values between 0
and 1.

Best regards,
Mike

raster.SGDF <- function (x, attr = 1, xcol = 1, ycol = 2, col = heat.colors(12),
    red = NULL, green = NULL, blue = NULL, axes = FALSE, xlim = NULL,
    ylim = NULL, add = FALSE, ..., asp = NA, setParUsrBB = FALSE,
interpolate = FALSE, angle = 0)
{
    ## function to scale values to [0, 1]
    scl <- function(x) (x - min(x, na.rm  = TRUE))/diff(range(x, na.rm = TRUE))
    ## bounding box of the image
    bb <- bbox(x)

    ## set up the background plot if it's not already
    if (!add)
        plot(as(x, "Spatial"), xlim = xlim, ylim = ylim, axes = axes,
            asp = asp, ..., setParUsrBB = setParUsrBB)

    ## 1-band case
    if (is.null(red)) {
        x <- x[attr]
        NAs <- is.na(x[[1]])
        nvalues <- length(unique(x[[1]][!NAs]))

        m <-  scl(t(matrix(x[[1]], x at grid@cells.dim[1], x at grid@cells.dim[2])))
        m <- matrix(col[as.vector(m) * (length(col)-1) + 1], nrow(m), ncol(m))
        ## if missing, set to white
        m[is.na(m)] <- rgb(1, 1, 1)


    } else {
        ## 3-band RGB case

        if (is.null(green) || is.null(blue))
            stop("all colour bands must be given")
        ## band data and missing values
        xd <- x at data[, c(red, green, blue)]
        NAs <- is.na(xd[, 1]) | is.na(xd[, 2]) | is.na(xd[, 3])
        if (any(NAs))
            xd <- xd[!NAs, ]
        ## create RGBs (using alpha=1 by default)
        RGBs <- rgb(xd, max = 255)
        if (any(NAs)) {
            z <- rep(NA, length(NAs))
            z[!NAs] <- RGBs
            RGBs <- z
        }

        cv <- coordinatevalues(getGridTopology(x))
        m <- t(matrix(RGBs, x at grid@cells.dim[1], x at grid@cells.dim[2],
        byrow = FALSE))
    }
    raster(m, bb[1,1], bb[2,1], bb[1,2], bb[2,2], interpolate =
interpolate, angle = angle)
}
op <- par(mfrow = c(3, 1))
raster.SGDF(Rlogo, red = "band1", green = "band1", blue = "band3",
interpolate = FALSE)

raster.SGDF(Rlogo, col = grey(seq(0, 1, length = 10)), angle = 28)
raster.SGDF(Rlogo, col = grey(seq(0, 1, length = 10)), interpolate = TRUE)
par(op)

 sessionInfo()
R version 2.11.0 Under development (unstable) (2010-03-07 r51225)
x86_64-pc-mingw64

locale:
[1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252
[3] LC_MONETARY=English_Australia.1252 LC_NUMERIC=C
[5] LC_TIME=English_Australia.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sp_0.9-60

loaded via a namespace (and not attached):
[1] grid_2.11.0    lattice_0.18-3 tools_2.11.0