From 2013a07b8c3fc07b5c7f99a532816363df464cc4 Mon Sep 17 00:00:00 2001 From: edzer Date: Mon, 6 May 2024 20:21:18 +0200 Subject: [PATCH] add gallery as vignette; #146 --- vignettes/csdacm.Rnw | 5 +- vignettes/sp_gallery.Rmd | 457 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 460 insertions(+), 2 deletions(-) create mode 100644 vignettes/sp_gallery.Rmd diff --git a/vignettes/csdacm.Rnw b/vignettes/csdacm.Rnw index d1f60ef..fef08f3 100644 --- a/vignettes/csdacm.Rnw +++ b/vignettes/csdacm.Rnw @@ -898,7 +898,7 @@ spplot.stgdf <- function(obj, zcol = 1, ..., format = NULL) { if (is.null(format)) format <- "%Y-%m-%d %H:%M:%S" cc <- coordinates(obj) df <- unstack(data.frame(obj[[zcol]], cc[,3])) - ns <- as.character(coordinatevalues(getGridTopology(obj))[[3]] + ISOdate(1970,1,1,0,0,0), format = format) + ns <- format(coordinatevalues(getGridTopology(obj))[[3]] + ISOdate(1970,1,1,0,0,0), format = format) cc2d <- cc[cc[,3] == min(cc[,3]), 1:2] obj <- SpatialPixelsDataFrame(cc2d, df) spplot(obj, names.attr = ns,...) @@ -913,7 +913,8 @@ setMethod("spplot", "SpatialTimeGridDataFrame", spplot.stgdf) .iwidth <- 6 .iheight <- 4 <> -print(spplot(x, format = "%H:%M:%S", as.table=TRUE)) +#print(spplot(x, format = "%H:%M:%S", as.table=TRUE)) +print(spplot(x, as.table=TRUE)) <> <> @ diff --git a/vignettes/sp_gallery.Rmd b/vignettes/sp_gallery.Rmd new file mode 100644 index 0000000..3fb3aff --- /dev/null +++ b/vignettes/sp_gallery.Rmd @@ -0,0 +1,457 @@ +--- +title: "sp Gallery" +output: + html_document: + toc: true + toc_float: + collapsed: false + smooth_scroll: false + toc_depth: 2 +vignette: > + %\VignetteIndexEntry{sp map gallery} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +# Plotting maps with `sp` + +This document shows example images created with objects +represented by one of the classes for spatial data in +packages sp. + +## 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 `sf`, and is loaded by + +```{r} +library(sf) +nc <- as(read_sf(system.file("gpkg/nc.gpkg", package="sf")), "Spatial") +``` + +## 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) + as(st_as_sf(m), "Spatial") +} +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)) + +``` + +## maps using ggplot2 + +(it is recommended to migrate to `sf`, and use `geom_sf()` for this) + +## 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",eval=FALSE} +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() +```