Chapter 9 Base and grid plots

With base plots we mean the plot methods as offered by base R. Higher-level plots created with ggplot2, which are built on top of the grid package, are discussed in chapter 10.

9.1 Base plots

The nice thing about base plots is that they work relatively simply, and can be built incrementally: after plotting a first element, it is possible to add further elements — like titles, axes, legends, annotations — with additional function calls. Making more complex plots such as faceted plots is more of a challenge with base plot. They may also be faster than e.g. ggplot2. This chapter focuses on the base plot methods for sf and stars objects, and in particular where they deviate from what one expects from base R plots: default colors, color key placement, and multiple maps.

9.1.1 plot.sf defaults

When we plot an sf object with multiple attributes, we get to see multiple attributes, as shown in figure 9.1; the warning indicates that the number of attributes plotted has been limited (by the argument max.plot), that by default no key is shown when multiple attributes are plotted, and that factor or character columns are plotted with a categorical scale, and numeric variables by a continuous color scale.

system.file("gpkg/nc.gpkg", package="sf") %>% read_sf -> nc
#> Warning: plotting the first 9 out of 14 attributes; use max.plot = 14 to plot
#> all
default plot with multiple attributes

Figure 9.1: default plot with multiple attributes

When the attributes values have a common reference system it makes sense to plot a key, and this can be done by

nc %>% select(SID74, SID79) %>% plot(key.pos = 4) # 1: below; 4: right
two attributes sharing a key

Figure 9.2: two attributes sharing a key

9.1.2 Controlling color, color breaks, key

The color of each feature geometry, be it a point, line or polygon, can be controlled by passing a single value to argument col, or to pass a vector with colors of length equal to the number of features. In that case, one does not get an automatic key, since the mapping of values to the key is unclear.

One can get further control over the key by specifying either of

  • breaks, to one of the classInt::classIntervals styles (“fixed”, “sd”, “equal”, “pretty”, “quantile”, “kmeans”, “hclust”, “bclust”, “fisher”, or “jenks”), or to the numeric vector of (increasing) class break values,
  • nbreaks to set the number of color breaks, and
  • pal, to pass a function that generates a palette, e.g. rainbow or viridis::viridis.
  • logz which causes legend values to be plotted as 10-powers
  • at which controls the values plotted along the key
  • key.pos plots the key beneath (1), left (2), over (3) or to the right (4) of the map(s), and is omitted if NULL.

An example is shown in figure 9.3.

plot(nc["SID74"], logz = TRUE, pal = viridis::viridis,
    breaks = c(0,.5,1,1.5,2,2.5), at = c(0,.5,1,1.5,2,2.5),
    key.width = lcm(1.3), key.length = 1)
controlling key color and breaks; log transform causes zero values to remain uncolored; custom key size parameters

Figure 9.3: controlling key color and breaks; log transform causes zero values to remain uncolored; custom key size parameters

For factor variables, the key shows factor levels, and key.length and/or key.width may need to be further controlled to get them look good.

9.1.3 Incrementally adding plot elements

We can add elements to maps without keys by simply adding add = TRUE. When a map has a key, the initial plot command needs to have reset = FALSE.

The reason for this is as follows. When plotting a map with a key on the side, in base plot one needs to cut up the plotting region in a region for the map, and a region for the key, similar to using par(mfrow = c(1,2)). In order to not keep the plot region splitting active for subsequent, unrelated plots, it is removed before the plot function returns. This means that by default one cannot add elements to a plot.

In case we want to add elements to a plot, we need to instruct plot not to reset, and then can add elements if we set add = TRUE, as shown in figure 9.4.

plot(nc["SID74"], pal = viridis::viridis, reset = FALSE, graticule = TRUE, axes = TRUE)
plot(st_point(c(-81.498,36.434)), col = 'red', pch = 1, cex = 4, lwd = 3, add = TRUE)
layout(matrix(1)) # manually reset the plotting region split-up
adding plot elements to a map with a key; graticule and axes

Figure 9.4: adding plot elements to a map with a key; graticule and axes

9.1.4 plotting graticules and axes

A graticule (section 8.5) is added if graticule = TRUE is set; figure 9.4 gives an example. Axes are by default omitted: they take space, and often other map elements, such as boundaries or coast lines, are sufficient for orientation. Axes can be added with axes = TRUE. They are given in values of the coordinate reference system of the data, unless a graticule is added, in which case they correspond to the graticule values.

9.1.5 plot.stars defaults

The base plot method for stars objects follows a number of parameters that were also discussed above for plot.sf. These are: key.pos, logz, axes, reset, key.width, key.length nbreaks, and breaks.

The default for breaks is "quantile"; along with the default join_zlim = TRUE this results in quantile values taken from all of the images shown in a composite plot to yeild color breaks. The result of this is shown in figure 4.1, the default stars plot for a multi-band raster. If join_zlim is FALSE, color breaks are computed separately for each subplot, and a joint legend is not possible, hence, no legend is drawn. For large datasets with more than 10000 pixels, a (regular) sample of size 10000 is taken to determine quantiles, rather than using the entire dataset. The col argument provides the pallette used, and not a color vector for each pixel; by default a linear grey scale is used. box_col is the color used for the boxes around subplots.

With the default of downsample = TRUE, plot.stars does not plot many more pixels than can be seen on the plotting device.

As shown in the left subfigure of 4.2, the rgb parameter can be used to specify three bands for creating rgb colors. maxColorValue, by default the maximum of all bands, can be used to specify the maximum pixel values for scaling pixel values into rgb components.

The effect of setting text_values and interpolate are shown in figure 9.5 showing a 10 x 10 cell raster. text_values prints cell values as text, e.g. for numeric examples. interpolate interpolates rgb values, and seems to limit its extent to the raster cell centres, rather than corners.

tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
par(mfrow = c(1, 2))
plot(x[,1:10,1:10,1], text_values=TRUE, key.pos = NULL, col = viridis::viridis(10), 
    reset = FALSE)
plot(x[,1:10,1:10,1:3], rgb=1:3, interpolate = TRUE, reset = FALSE)
stars plots with text values (left) and interpolated rgb colors (right); both maps plot a 10 x 10 raster

Figure 9.5: stars plots with text values (left) and interpolated rgb colors (right); both maps plot a 10 x 10 raster

Some rasters cannot be plotted by R’s image, but have to be plotted as features (polygons or points); examples are rotated, sheared and curvilinear grids. Figure 9.6 shows a sheared grid, and illustrates st_as_sf for conversion to simple features, merge for whether to merge polygons with identical values, and as_points to plot raster cells as points.

par(mfrow = c(1, 3))
xs = adrop(x[,1:10,1:10,1])
attr(attr(xs, "dimensions"), "raster")$affine = c(1, 3)
plot(xs, col = viridis::viridis(10), key.pos = NULL, reset = FALSE)
plot(st_as_sf(xs, as_points = FALSE, merge = FALSE), pal = viridis::viridis, 
    key.pos = NULL, reset = FALSE)
plot(st_as_sf(xs, as_points = TRUE, merge = FALSE), key.pos = NULL, 
    pal = viridis::viridis, pch = 16, reset = FALSE)
sheared rasters, plotted as simple features; left: with merged polygon features, middle: with seperate polygons, right: as point features

Figure 9.6: sheared rasters, plotted as simple features; left: with merged polygon features, middle: with seperate polygons, right: as point features

9.2 Combining base plots

Figure 9.6 is an exmaple of multiple plots of stars and sf objects in one plot. They could be combined because they lack the automatic legend key.

The explanation about the need for parameter reset and the examples above with multiple subplots setting it to FALSE already show the limitations of this approach: if we want a key next to each map, we cannot easily combine maps in subplots, using sf‘s or starsplot methods. If we can do without the automatic legend, all is fine.

A more flexible approach to combine maps and map elements into a composite is using package grid as shown in the next section, or ggplot facet plots (which also use grid) as will be shown in chapter 10 about ggplot2.

9.3 Grid plots and viewports

Package grid, one of the R base packages, takes a more structured approach to building plots than base plot does. It allows the creation of graphic objects (objects of class grob, or grobs) that contain all the information for plotting, and the definition of viewports, plotting subregions within which coordinate systems can easily be redefined. Packages like ggplot2 (chapter 10) and lattice are built upon the logic provided by grid.

Package sf contains st_as_grob methods for all feature geometry classes. Figure 9.7 shows a simple plot of all the nc counties drawn using grid functions.

system.file("gpkg/nc.gpkg", package="sf") %>% read_sf -> nc
st_viewport(nc) %>% pushViewport
st_geometry(nc) %>% lapply(st_as_grob) %>% lapply(grid.draw) -> x
simple plot created with grid functions

Figure 9.7: simple plot created with grid functions

The st_as_grob methods are exported by sf and used by ggplot2::geom_sf to convert geometries into plot-able objects; this way, ggplot2 needs no knowledge on how geometries are stored in R objects.