I have a raster file 'airtemp' and a polygon shapefile 'continents'. I'd like to superimpose the 'continents' on 'airtemp', so the boundary of 'continents' is visible on top of 'airtemp'. I plot the raster file by levelplot (lattice). I read the polygon by readShapeSpatial (maptools) first and then plot.
The problem is levelplot and plot have different scales. Plot tends to have smaller frame. Sorry I don't have a reproducible sample, but I feel this is a rather common issue for geophysicists. I've found a similar question here:
http://r.789695.n4.nabble.com/overlaying-a-levelplot-on-a-map-plot-td2019419.html
but I don't quite understand the solution.
Answer
You can overlay the shapefile using the +.trellis and layer
functions from the latticeExtra package (which is automatically
loaded with rasterVis).
library(raster)
library(rasterVis)
Let's build some data to play. You can skip this part if you already have a raster file and a shapefile.
library(maps)
library(mapdata)
library(maptools)
## raster
myRaster <- raster(xmn=-100, xmx=100, ymn=-60, ymx=60)
myRaster <- init(myRaster, runif)
## polygon shapefile
ext <- as.vector(extent(myRaster))
boundaries <- map('worldHires', fill=TRUE,
xlim=ext[1:2], ylim=ext[3:4],
plot=FALSE)
## read the map2SpatialPolygons help page for details
IDs <- sapply(strsplit(boundaries$names, ":"), function(x) x[1])
bPols <- map2SpatialPolygons(boundaries, IDs=IDs,
proj4string=CRS(projection(myRaster)))
Now you plot the raster file with rasterVis::levelplot, the
shapefile with sp::sp.polygons, and the overall graphic is produced
with +.trellis and layer.
levelplot(myRaster) + layer(sp.polygons(bPols))
sp.polygons uses a transparent color as default for fill, but you can change it:
levelplot(myRaster) + layer(sp.polygons(bPols, fill='white', alpha=0.3))
