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

[Oscar Perpiñán Lamigueiro]

Hi,

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

Try the example below.

Best,

Oscar. 

## vectors of latitude and longitude with non-NA values
longitude <- sample(seq(-130, -65, .1), 2000, replace=TRUE)
latitude <- sample(seq(25, 50, .1), 2000, replace=TRUE)
xmn <- min(longitude);
xmx <- max(longitude)
ymn <- min(latitude);
ymx <- max(latitude)

## Build a raster fill with NA
myrast <- raster(xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx, ncols=162, nrows=162)
## Extract the cell numbers corresponding to non-NA values
mycells <- cellFromXY(myrast, cbind(longitude, latitude))
## and fill them with the appropiate values
myrast[mycells] <- runif(length(longitude))

## Ready to plot...
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)
Waichler, Scott R writes:

> 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


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