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