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.perpi...@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. _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo