<!--
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{sp map gallery}
-->

# Plotting maps with `sp`

This  document shows example images created with objects
represented by one of the classes for spatial data in
packages sp. The [sp_gallery.Rmd](sp_gallery.Rmd) and
[sp_gallery.R](sp_gallery.R) files are available here, or on
[github](https://github.com/edzer/sp/tree/gh-pages). This html file
was generated from the R markdown file (sp_gallery.Rmd) by

```
library(rmarkdown)
render("sp_gallery.Rmd")
```

and the corresponding .R file by

```
library(knitr)
purl("sp_gallery.Rmd")
```

In case you load the .Rmd file in rstudio and press `Knit HTML`,
or alternatively load the .R file and run `Source with Echo`
(ctrl-shift-enter), you may see slightly different plots, which is
due to the different defaults for plotting device and dimension.

## Loading data and packages

The Meuse data set is loaded using a demo script in package `sp`,
```{r}
library(sp)
demo(meuse, ask = FALSE, echo = FALSE) # loads the meuse data sets
class(meuse)
```

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

```{r}
library(maptools)
library(rgdal) # for readOGR
nc <- readOGR(system.file("shapes/", package="maptools"), "sids")
proj4string(nc) <- CRS("+proj=longlat +datum=NAD27")
```

## Using base plot

The basic `plot` command on a `Spatial` object gives just the geometry, without axes:

```{r}
plot(meuse)
```

axes can be added when `axes=TRUE` is given; plot elements can be added incrementally:

```{r}
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))
```

For local projection systems, such as in this case

```{r}
proj4string(meuse)
```

the aspect ratio is set to 1, to make sure that one unit north
equals one unit east.  Even if the data are in long/lat degrees,
an aspect ratio is chosen such that in the middle of the plot one
unit north approximates one unit east. For small regions, this works
pretty well; also note the degree symbols in the axes tic marks:

```{r, fig.height=8}
crs.longlat = CRS("+init=epsg:4326")
meuse.longlat = spTransform(meuse, crs.longlat)
plot(meuse.longlat, axes = TRUE)
```

which looks different from the plot where one degree north (latitude) equals one
degree east (longitude):

```{r}
par(mar = rep(0,4))
plot(meuse.longlat, asp = 1)
```

## Graticules 

Instead of axes with ticks and tick marks, maps often have
graticules, a grid with constant longitude and latitude lines.
`sp` provides several helper functions to add graticules, either
in the local reference system, or in long/lat. Here is an
example of the local reference system:

```{r, fig.height=7.5}
par(mar = c(0, 0, 1, 0))
library(methods) # as
plot(as(meuse, "Spatial"), expandBB = c(.05, 0, 0, 0))
plot(gridlines(meuse), add = TRUE, col = grey(.8))
plot(meuse, add = TRUE)
text(labels(gridlines(meuse)), col = grey(.7))
title("default gridlines with Meuse projected data")
```

```{r}
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.
For 

```{r}
# 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))
pts=SpatialPoints(rbind(c(-180,-70),c(0,-70),c(180,-89),c(180,-70)), CRS("+init=epsg:4326"))
gl = gridlines(pts, easts = seq(-180,180,20), ndiscr = 100)
polar = CRS("+init=epsg:3031")
plot(spTransform(pts, polar), expandBB = c(.05, 0, .05, 0))
gl.polar = spTransform(gl, polar)
lines(gl.polar)
l = labels(gl.polar, crs.longlat, side = 3)
l$pos = NULL # pos is too simple, use adj:
text(l, adj = c(0.5, -0.5), cex = .8) 
l = labels(gl.polar, crs.longlat, side = 4)
l$srt = 0 # otherwise they end up upside-down
text(l, cex = .8)
title("grid line labels on polar projection, epsg 3031")

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")

```

## Other spatial objects in base plot

The following plot shows polygons with county name as labels at
their center point:

```{r}
par(mar = c(0, 0, 1, 0))
plot(nc)
invisible(text(coordinates(nc), labels=as.character(nc$NAME), cex=0.4))
```

This plot of a `SpatialPolygonsDataFrame` uses grey shades:

```{r}
names(nc)
rrt <- nc$SID74/nc$BIR74
brks <- quantile(rrt, seq(0,1,1/7))
cols <- grey((length(brks):2)/length(brks))
dens <- (2:length(brks))*3
par(mar = c(0, 0, 1, 0))
plot(nc, col=cols[findInterval(rrt, brks, all.inside=TRUE)])
```

The following plot shows a `SpatialPolygonsDataFrame`, using line densities

```{r}
rrt <- nc$SID74/nc$BIR74
brks <- quantile(rrt, seq(0,1,1/7))
cols <- grey((length(brks):2)/length(brks))
dens <- (2:length(brks))*3
par(mar = rep(0,4))
plot(nc, density=dens[findInterval(rrt, brks, all.inside=TRUE)])
```

Plot/image of a grid file, using base plot methods:

```{r}
image(meuse.grid)
```

```{r}
plot(meuse.grid["dist"])
points(meuse, col = 'green')
```

```{r}
plot(meuse.grid["dist"], zlim = c(0,1))
plot(geometry(meuse.grid), add = TRUE, col = grey(.8))
```

Read
[this](http://r-spatial.org/r/2016/03/08/plotting-spatial-grids.html)
blog post to find out more about the options available (and
limitations) for plotting gridded data with base plot methods in sp.

## using lattice plot (spplot)

The following plot colours points with a legend in the plotting area and adds scales:

```{r}
spplot(meuse, "zinc", do.log = TRUE,
	key.space=list(x = 0.1, y = 0.95, corner = c(0, 1)),
	scales=list(draw = TRUE))
```

The following plot has coloured points plot with legend in plotting area and scales;
it has a non-default number of cuts with user-supplied legend entries:
```{r}
spplot(meuse, "zinc", do.log = TRUE,
	key.space=list(x=0.2,y=0.9,corner=c(0,1)),
	scales=list(draw = TRUE), cuts = 3,
	legendEntries = c("low", "intermediate", "high"))
```

The following plot adds a scale bar and north arrow:

```{r}
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
	offset = c(178600,332490), scale = 500, fill=c("transparent","black"))
text1 = list("sp.text", c(178600,332590), "0")
text2 = list("sp.text", c(179100,332590), "500 m")
arrow = list("SpatialPolygonsRescale", layout.north.arrow(), 
	offset = c(178750,332000), scale = 400)
spplot(meuse, "zinc", do.log=T,
	key.space=list(x=0.1,y=0.93,corner=c(0,1)),
	sp.layout=list(scale,text1,text2,arrow),
	main = "Zinc (top soil)")
```

The following plot has north arrow and text outside panels

```{r}
rv = list("sp.polygons", meuse.riv, fill = "lightblue")
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
	offset = c(180500,329800), scale = 500, fill=c("transparent","black"), which = 1)
text1 = list("sp.text", c(180500,329900), "0", which = 1)
text2 = list("sp.text", c(181000,329900), "500 m", which = 1)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(), 
	offset = c(178750,332500), scale = 400)
spplot(meuse["zinc"], do.log = TRUE,
	key.space = "bottom", 
	sp.layout = list(rv, scale, text1, text2),
	main = "Zinc (top soil)",
	legend = list(right = list(fun = mapLegendGrob(layout.north.arrow()))))
```

The same plot; north arrow now inside panel, with custom panel function instead of sp.layout

```{r}
spplot(meuse, "zinc", panel = function(x, y, ...) {
		sp.polygons(meuse.riv, fill = "lightblue")
		SpatialPolygonsRescale(layout.scale.bar(), offset = c(179900,329600), 
			scale = 500, fill=c("transparent","black"))
		sp.text(c(179900,329700), "0")
		sp.text(c(180400,329700), "500 m")
		SpatialPolygonsRescale(layout.north.arrow(), 
			offset = c(178750,332500), scale = 400)
		panel.pointsplot(x, y, ...)
	},
	do.log = TRUE, cuts = 7,
	key.space = list(x = 0.1, y = 0.93, corner = c(0,1)),
	main = "Top soil zinc concentration (ppm)")
```

A multi-panel plot, scales + north arrow only in last plot: using
the `which` argument in a layout component (if `which=4` was set
as list component of sp.layout, the river would as well be drawn
only in that (last) panel)

```{r}
rv = list("sp.polygons", meuse.riv, fill = "lightblue")
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
	offset = c(180500,329800), scale = 500, fill=c("transparent","black"), which = 4)
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)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(), 
	offset = c(181300,329800), 
	scale = 400, which = 4)
cuts = c(.2,.5,1,2,5,10,20,50,100,200,500,1000,2000)
spplot(meuse, c("cadmium", "copper", "lead", "zinc"), do.log = TRUE,
	key.space = "right", as.table = TRUE,
	sp.layout=list(rv, scale, text1, text2, arrow), # note that rv is up front!
	main = "Heavy metals (top soil), ppm", cex = .7, cuts = cuts)
```

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

```{r}
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)
library(gstat)
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)
)
```

Reuse these results; universal kriging standard errors; grid plot
with point locations and polygon (river):

```{r}
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")
```

```{r}
arrow = list("SpatialPolygonsRescale", layout.north.arrow(),
    offset = c(-76,34), scale = 0.5, which = 2)
spplot(nc, c("SID74", "SID79"), names.attr = c("1974","1979"),
    colorkey=list(space="bottom"), scales = list(draw = TRUE),
    main = "SIDS (sudden infant death syndrome) in North Carolina",
    sp.layout = list(arrow), as.table = TRUE)
```

```{r}
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)
```

Bubble plots for cadmium and zinc:

```{r}
b1 = bubble(meuse, "cadmium", maxsize = 1.5, main = "cadmium concentrations (ppm)",
	key.entries = 2^(-1:4))
b2 = bubble(meuse, "zinc", maxsize = 1.5, main = "zinc concentrations (ppm)",
	key.entries =  100 * 2^(0:4))
print(b1, split = c(1,1,2,1), more = TRUE)
print(b2, split = c(2,1,2,1), more = FALSE)
```

Factor variables using `spplot`:

```{r}
# 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])
library(RColorBrewer)
## 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

```{r}
library(ggmap)
merc = CRS("+init=epsg:3857")
WGS84 = CRS("+init=epsg:4326")
meuse.ll = spTransform(meuse, WGS84)
bgMap = get_map(as.vector(bbox(meuse.ll)), source = "google", zoom = 13) # useless without zoom level
par(mar = rep(0,4))
plot(spTransform(meuse, merc), bgMap = bgMap, pch = 16, cex = .5)
```

Using spplot with ggmap-google bg:

```{r}
spplot(spTransform(meuse, merc), c("zinc",  "lead"), colorkey = TRUE,
	sp.layout = list(panel.ggmap, bgMap, first = TRUE))
```

Using open streetmap background from `ggmap`:

```{r}
par(mar = rep(0,4))
bb = t(apply(bbox(meuse.ll), 1, bbexpand, .04))
bgMap = get_map(as.vector(bb), source = "osm") # WGS84 for background map
plot(spTransform(meuse, merc), bgMap = bgMap, pch = 16, cex = .5)
```

Using `RgoogleMaps`:

```{r}
center = apply(coordinates(meuse.ll), 2, mean)[2:1]
library(RgoogleMaps)
g = GetMap(center=center, zoom=13) # google
par(mar = rep(0,4)) # fill full device
plot(spTransform(meuse, merc), bgMap = g, pch = 16, cex = .5)
```

Using `RgoogleMaps` and `spplot`:

```{r}
spplot(spTransform(meuse, merc), c("zinc",  "lead"), colorkey = TRUE,
	sp.layout = list(panel.RgoogleMaps, g, first = TRUE),
	scales = list(draw = TRUE))
```

Here is a Norway boundary example:

```{r}
# Norway boundary example:
library(cshapes)
cshp = cshp(as.Date("2000-01-1"))
norway = cshp[cshp$ISO1AL2 == "NO",]
NorwayMap = get_map(as.vector(bbox(norway))) # no is already in WGS84
par(mar = rep(0,4))
plot(spTransform(norway, merc), bgMap = NorwayMap, border = 'red')
```

## maps using ggplot2 

In order to plot polygon data with `ggplot` we need to prepare
(`fortify`) the data to a `data.frame`, and merge to each polygon
vertex (row) all corresponding attributes:

```{r}
library(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
```

next, we can plot the data, but if we want one distance unit north
to approximately match one unit east, we additionally need to set
`coord_map()` or `coord_quickmap()` to instruct `ggplot` this
concerns longitude/latitude coordinates (not stating the datum).
`coord_map()` projects to mercator (by default), `coord_quickmap()`
computes a suitable aspect ratio and plots degrees.

```{r}
ggplot(ff) + 
	aes(long, lat, group = group, fill = SID74) + 
	geom_polygon() + 
	coord_quickmap()
```

Bubble plots can be generated through `ggplot` by

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

where `coord_equal()` sets the unit aspect ratio of the axes,
for projected coordinates.

The [ggplot2 wiki](https://github.com/hadley/ggplot2/wiki/plotting-polygon-shapefiles)
indicates that plotting polygons with holes may be problematic with ggplot.

## Other packages for making maps

The [spatial task
view](https://cran.r-project.org/web/views/Spatial.html) contains
pointers to other CRAN packages that have functionality for creating
general or custom maps, in addition to pointers to packages that
can be used for analyzing spatial data.

## Interactive maps: leaflet, mapview

R packages leaflet and mapview provide interactive, browser-based
maps building upon the leaflet javascript library. Example with
points, grid and polygons follow:

```{r, results="markup"}
library(mapview)
mapview(meuse, zcol = c("zinc", "lead"), legend = TRUE)
```

```{r, results="markup"}
mapview(meuse.grid, zcol = c("soil", "dist"), legend = TRUE)
```

```{r, results="markup"}
mapview(nc, zcol = c("SID74", "SID79"), alpha.regions = 1.0, legend = TRUE)
```

Mapview also allows grids of view that are synced
```{r, results="markup"}
m1 <- mapview(meuse, zcol = "soil", burst = TRUE, legend = TRUE)
m2 <- mapview(meuse, zcol = "lead", legend = TRUE)
m3 <- mapview(meuse, zcol = "landuse", map.types = "Esri.WorldImagery", legend = TRUE)
m4 <- mapview(meuse, zcol = "dist.m", legend = TRUE)
sync(m1, m2, m3, m4) # 4 panels synchronised
# latticeView(m1, m2, m3, m4) # 4 panels
```

more examples are found [here](https://environmentalinformatics-marburg.github.io/web-presentations/20150723_mapView.html).

## SessionInfo

```{r}
sessionInfo()
```
