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

[Michael Sumner]

My apologies, I left out the data set up for Rlogo. Here is the full session:


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)
}
  library(sp)  ## R 2.11.0
  data(Rlogo)
  d = dim(Rlogo)
  cellsize = abs(c(gt[2],gt[6]))
  cells.dim = c(d[1], d[2]) # c(d[2],d[1])
  cellcentre.offset = c(x = gt[1] + 0.5 * cellsize[1], y = gt[4] -
(d[2] - 0.5) * abs(cellsize[2]))
  grid = GridTopology(cellcentre.offset, cellsize, cells.dim)
  df = as.vector(Rlogo[,,1])
  for (band in 2:d[3]) df = cbind(df, as.vector(Rlogo[,,band]))
  df = as.data.frame(df)
  names(df) = paste("band", 1:d[3], sep="")
  Rlogo <- SpatialGridDataFrame(grid = grid, data = df)

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




On Tue, Mar 9, 2010 at 1:13 PM, Michael Sumner <mdsumner at gmail.com> wrote:
> 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
>