I want to plot the country borders of North America over a raster image depicting some variable and then overlay contours on top of the plot using R. I have been successful in doing this using base graphics and lattice, but it seems that the plotting process is much too slow! I have not done this in ggplot2 yet, but I doubt that it will fare better in terms of speed.
I have the data in a netcdf file created from a grib file. For now, I downloaded the country borders for Canada, USA and Mexico, which were available in RData files from GADM which read into R as SpatialPolygonsDataFrame objects.
Here is some code:
# Load packages
library(raster)
#library(ncdf) # If you cannot install ncdf4
library(ncdf4)
# Read in the file, get the 13th layer
# fn <- 'path_to_file'
r <- raster(fn, band=13)
# Set the projection and extent
p4 <- "+proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +x_0=32.46341 +y_0=32.46341 +lon_0=-107 +lat_0=1.0"
projection(r) <- CRS(p4)
extent(r) <- c(-5648.71, 5680.72, 1481.40, 10430.62)
# Get the country borders
# This will download the RData files to your working directory
can<-getData('GADM', country="CAN", level=1)
usa<-getData('GADM', country="USA", level=1)
mex<-getData('GADM', country="MEX", level=1)
# Project to model grid
can_p <- spTransform(can, CRS(p4))
usa_p <- spTransform(usa, CRS(p4))
mex_p <- spTransform(mex, CRS(p4))
### USING BASE GRAPHICS
par(mar=c(0,0,0,0))
# Plot the raster
bins <- 100
plot(r, axes=FALSE, box=FALSE, legend=FALSE,
col=rev( rainbow(bins,start=0,end=1) ),
breaks=seq(4500,6000,length.out=bins))
plot(r, legend.only=TRUE, col=rev( rainbow(bins,start=0,end=1)),
legend.width=0.5, legend.shrink=0.75,
breaks=seq(4500,6000,length.out=bins),
axis.args=list(at=seq(4500,6000,length.out=11),
labels=seq(4500,6000,length.out=11),
cex.axis=0.5),
legend.args=list(text='Height (m)', side=4, font=2,
line=2, cex=0.8))
# Plot the borders
# These are so slow!!
plot(can_p, add=TRUE, border='white', lwd=2)
plot(usa_p, add=TRUE, border='white', lwd=2)
plot(mex_p, add=TRUE, border='white', lwd=2)
# Add the contours
contour(r, add=TRUE, nlevel=5)
### USING LATTICE
library(rasterVis)
# Some settings for our themes
myTheme <- RdBuTheme()
myTheme$axis.line$col<-"transparent"
myTheme$add.line$alpha <- 1
myTheme2 <- myTheme
myTheme2$regions$col <- 'transparent'
myTheme2$add.text$cex <- 0.7
myTheme2$add.line$lwd <- 1
myTheme2$add.line$alpha <- 0.8
# Get JUST the contour lines
contours <- contourplot(r, margin=FALSE, scales=list(draw=FALSE),
par.settings=myTheme2, pretty=TRUE, key=NULL, cuts=5,
labels=TRUE)
# Plot the colour
levels <- levelplot(r, contour=FALSE, margin=FALSE, scales=list(draw=FALSE),
par.settings = myTheme, cuts=100)
# Plot!
levels +
layer(sp.polygons(can_p, col='green', lwd=2)) +
layer(sp.polygons(usa_p, col='green', lwd=2)) +
layer(sp.polygons(mex_p, col='green', lwd=2)) +
contours
Is there a way to speed up the plotting of the polygons?
On the system that I am working on, the plotting takes several minutes. I want to eventually make a function that will easily generate a number of these plots for inspection, and I reckon I will be plotting many of these maps, so I want to increase the speed of the plots!