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...@grid@cells.dim[1],
x...@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...@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...@grid@cells.dim[1], x...@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 <mdsum...@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...@grid@cells.dim[1],
> x...@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...@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...@grid@cells.dim[1], x...@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
>
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
