[R-sig-Geo] how to use levelplot() with a geographic projection

[Oscar Perpiñán Lamigueiro]

Hello,

Scott, here I propose a different way to build and fill the RasterStack
object:

## make the stack object
myrast <- raster(xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx, ncols=num.lon, nrows=num.lat)
myrast[] <- NA
raster.stack <- stack(lapply(1:4, function(i)myrast))

## use cellFromXY() to select the target cells in the raster
select.cell.numbers <- cellFromXY(raster.stack, select.coords)

## give the selected cells some values
vals <- sapply(1:4, function(i)seq(i, i+1, length=num.sel.cells))
raster.stack[select.cell.numbers] <- vals

Best,

Oscar

Waichler, Scott R writes:

> 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.


-- 
Oscar Perpiñán Lamigueiro
Grupo de Sistemas Fotovoltaicos (IES-UPM)
Dpto. Ingeniería Eléctrica (EUITI-UPM)
URL: http://http://oscarperpinan.github.io
Twitter: @oscarperpinan