Hi again,
In working on this some more I discovered that I need to transpose a data
matrix, then flip the raster object created from the matrix to get the proper
rendering in a plot. Can someone explain why that is?
My example has to do with the state of California, with 2918 1/8-degree cells.
I regret I can't provide a small working example--I don't know to make one up
that would clearly show correct vs. incorrect geographic rendering.
library(maps)
library(maptools)
library(raster)
library(rasterVis)
data(stateMapEnv)
state.map <- map("state", plot=F)
# longitude and latitude are vectors giving the respective coordinates
# for cells that are non-NA (n=2918)
xmn <- min(longitude); xmx <- max(longitude)
ymn <- min(latitude); ymx <- max(latitude)
x.bounds <- seq(xmn - 1/16, xmx + 1/16, by = 1/8) # grid spacing is 1/8 degree
y.bounds <- seq(ymn - 1/16, ymx + 1/16, by = 1/8)
ii <- findInterval(longitude, x.bounds) # indices
jj <- findInterval(latitude, y.bounds) # indices
index.matrix <- cbind(ii,jj)
num.long <- length(x.bounds) - 1
num.lat <- length(y.bounds) - 1
mat <- array(NA, dim=c(num.long, num.lat)) # full, enclosing array with
rows=long and columns=lat
mat[index.matrix] <- runif(n=length(longitude)) # insert non-NA values in the
right cells
# Why do I need to transpose the matrix then flip the raster for proper
rendering in levelplot?
myrast <- flip(raster(t(mat), xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx,
crs="+proj=longlat +datum=WGS84"),
direction="y")
ext <- as.vector(extent(myrast))
boundaries <- map('state', fill=TRUE, xlim=ext[1:2], ylim=ext[3:4], plot=FALSE)
IDs <- sapply(strsplit(boundaries$names, ":"), function(x) x[1])
bPols <- map2SpatialPolygons(boundaries, IDs=IDs,
proj4string=CRS(projection(myrast)))
levelplot(myrast,
panel = function(...) {
panel.levelplot(...)
sp.polygons(bPols)
},
margin=F)
--Scott Waichler
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
