how to use levelplot() with a geographic projection

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

-----Original Message-----
From: Oscar Perpi??n Lamigueiro [mailto:oscar.perpinan at gmail.com]
Subject: Re: [R-sig-Geo] how to use levelplot() with a geographic
projection

Hi,

I think your problem is related with the way you define the RasterLayer.
You should use cellFromXY instead.

For the archive, here is a complete example of using rasterVis to do the following:

-- use levelplot() to draw multiple panels with maps being geographically correct
-- create raster objects and supply values for some of the cells, selected from
   a matrix of longitude, latitude coordinates
-- overlay polygons (state boundaries) for reference

Thanks to Oscar Perpi??n Lamigueiro and Tom Philippi for their help.

library(maps)
library(maptools)
library(raster)
library(rasterVis)

data(stateMapEnv)  # for U.S. state boundaries
state.map <- map("state", plot=F)
set_Polypath(FALSE)  # per Roger Bivand, allows overprinting with sp.polygons()

# generate geographic cell boundaries that cover the northwest U.S.
lon.cell.bounds <- seq(-125, -109, by=1)
lat.cell.bounds <- seq(41, 49, by=1)
# save min and max for raster later
xmn <- min(lon.cell.bounds); xmx <- max(lon.cell.bounds)
ymn <- min(lat.cell.bounds); ymx <- max(lat.cell.bounds)
# get the centers of the cells
lon <- (lon.cell.bounds[1:(length(lon.cell.bounds)-1)] + lon.cell.bounds[2:length(lon.cell.bounds)]) / 2
lat <- (lat.cell.bounds[1:(length(lat.cell.bounds)-1)] + lat.cell.bounds[2:length(lat.cell.bounds)]) / 2
num.lon <- length(lon)  # number of cells in longitude (x)
num.lat <- length(lat)  # number of cells in latitude (y)

# Create a matrix of selected cell coordinates, roughly corresponding to the state of Washington.
select.lon <- seq(-123.5, -117.5, by=1)
select.lat <- seq(46.5, 48.5, by=1)
select.coords <- expand.grid(select.lon, select.lat)
num.sel.cells <- nrow(select.coords)

raster.layers <- list()
for(i in 1:4) {
  # make the raster object
  myrast <- raster(xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx, ncols=num.lon, nrows=num.lat)
  # use cellFromXY() to select the target cells in the raster
  select.cell.numbers <- cellFromXY(myrast, select.coords)
  # give the selected cells some values
  myrast[select.cell.numbers] <- seq(i, i+1, length=num.sel.cells)
  # save the raster to the list
  raster.layers[[i]] <- myrast  
}
raster.stack <- stack(raster.layers)  # stack is a collection of RasterLayer objects

# polygon shapefile for the state boundaries
ext <- as.vector(extent(raster.stack[[1]]))
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(raster.stack[[1]])))

# plot it!
levelplot(raster.stack,
          names.attr=c("Panel 1, 1-2", "Panel 2, 2-3", "Panel 3, 3-4", "Panel 4, 4-5"), 
          panel = function(...) {
              panel.levelplot(...)
              sp.polygons(bPols)
          }, main=paste(sep="\n", "Washington should be mostly covered with a rectangle",
               "and color should be darkest in lower left corner,",
               "lightest in upper right"),
          margin=F)

-----Original Message-----
From: Oscar Perpi??n Lamigueiro [mailto:oscar.perpinan at gmail.com]
Subject: Re: [R-sig-Geo] how to use levelplot() with a geographic
projection

Hi,

I think your problem is related with the way you define the RasterLayer.
You should use cellFromXY instead.