loadLibs = function(libs="RgoogleMaps"){
  for (l in libs) 
    suppressMessages(suppressWarnings(require(l,character.only =TRUE, quietly=TRUE)))
}
if (!require(RgoogleMaps))  install.packages("RgoogleMaps", repos="http://R-Forge.R-project.org")

pckgsForExercises = c("geosphere", "ggmap", "RgoogleMaps", "maps", "gstat", "ggplot2")
pckgs2Compile = c("sp", "rgdal", "leaflet","maptools","mapview", "rgeos", "RColorBrewer")

pckgs2load = c(pckgsForExercises,pckgs2Compile)

for (pkg in pckgs2load) loadLibs(pkg) 
#The Meuse data set is loaded using a demo script in package `sp`,
library(sp)
demo(meuse, ask = FALSE, echo = FALSE) # loads the meuse data sets
class(meuse)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
crs.longlat = CRS("+init=epsg:4326")
meuse.longlat = spTransform(meuse, crs.longlat)

#The North Carolina SIDS (sudden infant death syndrome) data set
#is available from package `maptools`, and is loaded by

nc <- readOGR(system.file("shapes/", package="maptools"), "sids")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:/user/home/loecherm/Documents/R/win-library/3.3/maptools/shapes", layer: "sids"
## with 100 features
## It has 14 fields
## Integer64 fields read as strings:  CNTY_ CNTY_ID FIPSNO
proj4string(nc) <- CRS("+proj=longlat +datum=NAD27")

OFFLINE=FALSE # use previously downloaded maps/data
set.seed(123)

base plot

proj4string(meuse)
## [1] "+init=epsg:28992 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs"
par(mar=c(7,2,1,1))
plot(meuse, pch = 1, cex = sqrt(meuse$zinc)/12, axes = TRUE)
v = c(100,200,400,800,1600)
legend("topleft", legend = v, pch = 1, pt.cex = sqrt(v)/12)
plot(meuse.riv, add = TRUE, col = grey(.9, alpha = .5))

the aspect ratio is set to 1, to make sure that one unit north equals one unit east.

Graticules

graticules: a grid with constant longitude and latitude lines.

par(mar = c(0, 0, 1, 0))
grd <- gridlines(meuse.longlat)
grd_x <- spTransform(grd, CRS(proj4string(meuse)))
plot(as(meuse, "Spatial"), expandBB = c(.05, 0, 0, 0))
plot(grd_x, add=TRUE, col = grey(.8))
plot(meuse, add = TRUE)
text(labels(grd_x, crs.longlat), col = grey(.7))
title("longitude latitude gridlines and labels")

These lines look pretty straight, because it concerns a small area.

Larger Area

# demonstrate axis labels with angle, both sides:

maps2sp = function(xlim, ylim, l.out = 100, clip = TRUE) {
    stopifnot(require(maps))
    m = map(xlim = xlim, ylim = ylim, plot = FALSE, fill = TRUE)
    p = rbind(cbind(xlim[1], seq(ylim[1],ylim[2],length.out = l.out)),
          cbind(seq(xlim[1],xlim[2],length.out = l.out),ylim[2]),
          cbind(xlim[2],seq(ylim[2],ylim[1],length.out = l.out)),
          cbind(seq(xlim[2],xlim[1],length.out = l.out),ylim[1]))
    LL = CRS("+init=epsg:4326")
    bb = SpatialPolygons(list(Polygons(list(Polygon(list(p))),"bb")), proj4string = LL)
    IDs <- sapply(strsplit(m$names, ":"), function(x) x[1])
    stopifnot(require(maptools))
    m <- map2SpatialPolygons(m, IDs=IDs, proj4string = LL)
    if (!clip)
        m
    else {
        stopifnot(require(rgeos))
        gIntersection(m, bb) # cut map slice in WGS84
    }
}
par(mar = c(0, 0, 1, 0))
m = maps2sp(c(-100,-20), c(10,55))
sp = SpatialPoints(rbind(c(-101,9), c(-101,55), c(-19,9), c(-19,55)), CRS("+init=epsg:4326"))
laea = CRS("+proj=laea +lat_0=30 +lon_0=-40")
m.laea = spTransform(m, laea)
sp.laea = spTransform(sp, laea)
plot(as(m.laea, "Spatial"), expandBB = c(.1, 0.05, .1, .1))
plot(m.laea, col = grey(.8), add = TRUE)
gl = gridlines(sp, easts = c(-100,-80,-60,-40,-20), norths = c(20,30,40,50))
gl.laea = spTransform(gl, laea)
plot(gl.laea, add = TRUE)
text(labels(gl.laea, crs.longlat))
text(labels(gl.laea, crs.longlat, side = 3:4), col = 'red')
title("curved text label demo")

polar:

par(mar = c(0, 0, 1, 0))
m = maps2sp(xlim = c(-180,180), ylim = c(-90,-70), clip = FALSE)
gl = gridlines(m, easts = seq(-180,180,20))
polar = CRS("+init=epsg:3031")
gl.polar = spTransform(gl, polar)
plot(as(gl.polar, "Spatial"), expandBB = c(.05, 0, .05, 0))
plot(gl.polar, add = TRUE)
plot(spTransform(m, polar), add = TRUE, col = grey(0.8, 0.8))
l = labels(gl.polar, crs.longlat, side = 3)
# pos is too simple here, use adj:
l$pos = NULL 
text(l, adj = c(0.5, -0.3), cex = .8)
l = labels(gl.polar, crs.longlat, side = 2)
l$srt = 0 # otherwise they are upside-down
text(l, cex = .8)
title("grid line labels on polar projection, epsg 3031")

Finally some geostatistics

Comparing four kriging varieties in a multi-panel plot with shared scale:

rv = list("sp.polygons", meuse.riv, fill = "blue", alpha = 0.1)
pts = list("sp.points", meuse, pch = 3, col = "grey", alpha = .5)
text1 = list("sp.text", c(180500,329900), "0", cex = .5, which = 4)
text2 = list("sp.text", c(181000,329900), "500 m", cex = .5, which = 4)
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
    offset = c(180500,329800), scale = 500, fill=c("transparent","black"), which = 4)

v.ok = variogram(log(zinc)~1, meuse)
ok.model = fit.variogram(v.ok, vgm(1, "Exp", 500, 1))
# plot(v.ok, ok.model, main = "ordinary kriging")
v.uk = variogram(log(zinc)~sqrt(dist), meuse)
uk.model = fit.variogram(v.uk, vgm(1, "Exp", 300, 1))
# plot(v.uk, uk.model, main = "universal kriging")
meuse[["ff"]] = factor(meuse[["ffreq"]])
meuse.grid[["ff"]] = factor(meuse.grid[["ffreq"]])
v.sk = variogram(log(zinc)~ff, meuse)
sk.model = fit.variogram(v.sk, vgm(1, "Exp", 300, 1))
# plot(v.sk, sk.model, main = "stratified kriging")
zn.ok = krige(log(zinc)~1,          meuse, meuse.grid, model = ok.model, debug.level = 0)
zn.uk = krige(log(zinc)~sqrt(dist), meuse, meuse.grid, model = uk.model, debug.level = 0)
zn.sk = krige(log(zinc)~ff,         meuse, meuse.grid, model = sk.model, debug.level = 0)
zn.id = krige(log(zinc)~1,          meuse, meuse.grid, debug.level = 0)

zn = zn.ok
zn[["a"]] = zn.ok[["var1.pred"]]
zn[["b"]] = zn.uk[["var1.pred"]]
zn[["c"]] = zn.sk[["var1.pred"]]
zn[["d"]] = zn.id[["var1.pred"]]
spplot(zn, c("a", "b", "c", "d"), 
    names.attr = c("ordinary kriging", "universal kriging with dist to river", 
        "stratified kriging with flood freq", "inverse distance"), 
    as.table = TRUE, main = "log-zinc interpolation",
    sp.layout = list(rv, scale, text1, text2)
)

universal kriging standard errors

rv = list("sp.polygons", meuse.riv, fill = "blue", alpha = 0.1)
pts = list("sp.points", meuse, pch = 3, col = "grey", alpha = .7)
# spplot(zn.uk, "var1.pred",
#   sp.layout = list(rv, scale, text1, text2, pts),
#   main = "log(zinc); universal kriging using sqrt(dist to Meuse)")
zn.uk[["se"]] = sqrt(zn.uk[["var1.var"]])
 spplot(zn.uk, "se", sp.layout = list(rv, pts),
    main = "log(zinc); universal kriging standard errors")

SID

arrow = list("SpatialPolygonsRescale", layout.north.arrow(), 
    offset = c(-76,34), scale = 0.5, which = 2)
#scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
#   offset = c(-77.5,34), scale = 1, fill=c("transparent","black"), which = 2)
#text1 = list("sp.text", c(-77.5,34.15), "0", which = 2)
#text2 = list("sp.text", c(-76.5,34.15), "1 degree", which = 2)
# create a fake lines data set:
## multi-panel plot with coloured lines: North Carolina SIDS
spplot(nc, c("SID74","SID79"), names.attr = c("1974","1979"),
    colorkey=list(space="bottom"),
    main = "SIDS (sudden infant death syndrome) in North Carolina",
    sp.layout = arrow, as.table = TRUE)

Factor variables using spplot:

# create two dummy factor variables, with equal labels:
set.seed(31)
nc$f = factor(sample(1:5, 100,replace = TRUE),labels=letters[1:5])
nc$g = factor(sample(1:5, 100,replace = TRUE),labels=letters[1:5])

## Two (dummy) factor variables shown with qualitative colour ramp; degrees in axes
spplot(nc, c("f","g"), col.regions=brewer.pal(5, "Set3"), scales=list(draw = TRUE))

Adding web map backgrounds

Both base plots and spplot can add web map backgrounds to plots, to help referencing. Special care needs to be taken to project the data into the web mercator (epsg 3857) projection.

Using regular plot and ggmap background maps

merc = CRS("+init=epsg:3857")
WGS84 = CRS("+init=epsg:4326")
meuse.ll = spTransform(meuse, WGS84)
bgMap = GetMap(...)

plot(spTransform(meuse, merc), bgMap = bgMap)

RgoogleMaps and spplot:

center = apply(coordinates(meuse.ll), 2, mean)[2:1]
g = GetMap(center=center, zoom=13) # google
spplot(spTransform(meuse, merc), c("zinc",  "lead"), colorkey = TRUE,
    sp.layout = list(panel.RgoogleMaps, g, first = TRUE),
    scales = list(draw = TRUE))

maps using ggplot2

f = fortify(nc, region = "CNTY_ID") 
ff = merge(f, nc@data, by.x = "id", by.y = "CNTY_ID") # this replicates a lot of info
ggplot(ff) + 
    aes(long, lat, group = group, fill = SID74) + 
    geom_polygon() + 
    coord_quickmap()

Bubble plots in ggplot

ggplot(as.data.frame(meuse)) + 
    aes(x, y, size = zinc) + 
    geom_point() +
    coord_equal()

Next: leaflet and mapview