Plotting vector maps overlaid on top of a raster background in R

Agustin.Lobo@ictja.csic.es
20140507
updated: 20140513
file: /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/vectorOnraster_log.R

Objects returned by dismo:gmap() are raster layers of package raster with the correct CRS information (“Google PseudoMercator” epsg:3857) that can be used as background for displaying other geographic layers in R provided their CRS are consistent. We show here the different alternatives.

Requirements

require(utils)
require(rgdal)
require(raster)
require(dismo)
require(rasterVis)

1. Input data

1.1 Vector data

download.file("https://dl.dropboxusercontent.com/u/3180464/NWEur.zip", "NWEur.zip", 
    method = "curl")
unzip("NWEur.zip", exdir = "NWEur")
NWEur <- readOGR(dsn = "./NWEur", layer = "NWEur", stringsAsFactors = FALSE)

## OGR data source with driver: ESRI Shapefile 
## Source: "./NWEur", layer: "NWEur"
## with 9 features and 18 fields
## Feature type: wkbPolygon with 2 dimensions

projection(NWEur)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

1.2 Background raster layer

migmap <- gmap(x = NWEur, type = "hybrid", zoom = 6, file = "NWEurdismo")
migmap

## class       : RasterLayer 
## dimensions  : 640, 378, 241920  (nrow, ncol, ncell)
## resolution  : 2446, 2446  (x, y)
## extent      : 945288, 1869870, 5982446, 7547876  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
## data source : /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/NWEurdismo.gif 
## names       : NWEurdismo 
## values      : 0, 255  (min, max)

projection(NWEur)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

1.3 Consistent CRS

As the CRS of each layer is different, we reproject NWEur to the CRS of the gmap object:

NWEurGM <- spTransform(NWEur, CRS(projection(migmap)))

2 Standard plots

Overlay using standard plots: note the vector layer must have the same CRS as the raster layer

plot(migmap)

plot(NWEurGM, add = TRUE, lwd = 2, ext = extent(migmap), border = "yellow")

Note you must decrease the zoom level in gmap() in case the extent is too large: google does not let you download large extents at high resolution. The following is closer to the extent of NWEur:

migmap <- gmap(x = NWEur, type = "hybrid", zoom = 4, file = "NWEurdismo")
plot(migmap)

plot(NWEurGM, add = TRUE, lwd = 2, ext = extent(migmap), border = "yellow")

3. Using rasterVis

3.1 Using levelplot() to plot the raster layer and layer() to overlay the reference polygons layer

colors <- migmap@legend@colortable
levelplot(migmap, maxpixels = ncell(migmap), col.regions = colors, at = 0:255, 
    panel = panel.levelplot.raster, interpolate = TRUE, colorkey = FALSE, margin = FALSE) + 
    layer(sp.polygons(NWEurGM, col = "yellow"))

3.2 Using levelplot() to plot the raster layer and spplot() to overlay the reference polygons layer

levelplot(migmap, maxpixels = ncell(migmap), col.regions = colors, at = 0:255, 
    panel = panel.levelplot.raster, interpolate = TRUE, colorkey = FALSE, margin = FALSE) + 
    spplot(NWEurGM, zcol = "NAME", fill = "transparent", col = "yellow", xlim = c(extent(migmap)@xmin, 
        extent(migmap)@xmax), ylim = c(extent(migmap)@ymin, extent(migmap)@ymax), 
        colorkey = FALSE)