1 Plotting in quickPlot

Most visualization approaches in R (e.g., base, RgoogleVis, ggplot2, rgl) are not optimized for speed or modularity. Rather, they are designed to build publication-quality or interactive figures that do not need to be redrawn hundreds or more times, maybe even in rapid succession. Plotting requirements for many other uses, including simulation modelling, are different. Changes to data, maps, objects can happen rapidly and repeatedly. Visualizing this can be critical, especially during the model building phase. Furthermore, no plotting method has been designed for modularity – the ability to plot something without being aware of what is already there.

quickPlot comes with a plotting function to fulfill these two gaps. The main plotting function, Plot (i.e., with a capital P), is built using the grid package. We have specifically built a plotting system that allows for relatively fast plotting of rasters (in the raster package), SpatialPoints, SpatialPolygons, SpatialLines (all within the sp package), ggplot objects, histogram, and igraph (package igraph) objects with the ability to make multi-frame plots without knowing which plots are already plotted. In other words, the main plotting function can handle modular plotting, where each call to Plot adds plots, without each knowing what the current state of the active plotting device is. This means that plotting can be treated as modular. The Plot function is therefore intended to be used as a way to interact visually especially during development of rapid, repeated visualizations, such as simulation modelling. If fine tuning and customization are desired, other plotting tools may be more suited (e.g., ggplot2, or a dedicated GIS program).

To demonstrate plotting, we first load some maps. These maps are randomly generated maps that come with the quickPlot package. In the code snippet below, we create the list of files to load, which is every file in the maps/ subdirectory of the package. Then we load that list of files.

#  Make list of maps from package database to load, and what functions to use to load them
library(data.table)
library(terra)
library(quickPlot)

# omit forestAge and percentPine maps for simplicity
files <- system.file("maps", package = "quickPlot") |>
  dir(full.names = TRUE, pattern = "tif")
filelist <- data.frame(file = files[-c(2, 5) ], stringsAsFactors = FALSE)
print(filelist)

##                                                                file
## 1            /tmp/Rtmpjnkjra/Rinst57c35d7140/quickPlot/maps/DEM.tif
## 2    /tmp/Rtmpjnkjra/Rinst57c35d7140/quickPlot/maps/forestCover.tif
## 3 /tmp/Rtmpjnkjra/Rinst57c35d7140/quickPlot/maps/habitatQuality.tif

# Load files to memory (using rasterToMemory), assign to a simList we call maps 
maps <- lapply(filelist$file, function(x) rast(x))
names(maps) <- sapply(basename(filelist$file), function(x) {
  strsplit(x, split = "\\.")[[1]][1]
})

# put into a single stack object in the simulation environment for ease of use below
landscape <- c(maps$DEM, maps$forestCover, maps$habitatQuality)

2 The Plot function

There are several features of Plot that are worth highlighting. A primary purpose of Plot is to plot spatial type objects; however, as noted below, other object types can also be plotted. Spatial objects can be layered on top of one another; non-spatial objects can not.

Plot(landscape, new = TRUE)

# make a SpatialPoints object
caribou <- terra::vect(cbind(x = stats::runif(1e2, -50, 50),
                             y = stats::runif(1e2, -50, 50)))
Plot(caribou)
Plot(caribou, addTo = "landscape$habitatQuality")

# from ?vect help file
x1 <- rbind(c(-18,-2), c(-14,5.5), c(1, 0), c(-14,-6))
x2 <- rbind(c(-1,0), c(14,6), c(16,0), c(14,-5.5))
x3 <- rbind(c(-12.5,0), c(0,6), c(4,0.5), c(1.5,-4.5))
hole <- rbind(c(8,0), c(10.5,1.3), c(12,0.2), c(10.5,-1.3))
z <- rbind(cbind(object=1, part=1, x1, hole=0), cbind(object=2, part=1, x3, hole=0),
            cbind(object=3, part=1, x2, hole=0), cbind(object=3, part=1, hole, hole=1))
colnames(z)[3:4] <- c('x', 'y')
spP <- terra::vect(z, "polygons")

Plot(spP)
Plot(spP, addTo = "landscape$habitatQuality", gp = gpar(lwd = 2))



# from ?vect help file
z[z[, "hole"]==1, "object"] <- 4
sl <- vect(z[,1:4], "lines")

Plot(sl, gp = gpar(col = c("red", "blue"), lwd = 2), addTo = "landscape$DEM")

clearPlot()
Plot(landscape, caribou, maps$DEM, spP, axes = TRUE,
     gp = gpar(cex = 0.5), visualSqueeze = 0.65)

if (requireNamespace("ggplot2")) {
  ggObj <- data.frame(x = stats::rnorm(1e3)) |> ggplot2::ggplot(ggplot2::aes(x = x)) +
    ggplot2::geom_histogram()
  clearPlot()
  Plot(caribou, axes = "L", new = TRUE) 
  Plot(ggObj) 
}

## Loading required namespace: ggplot2

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

baseObj <- rnorm(1e3)
baseObj2 <- baseObj * 1.2 + rnorm(1e3)
clearPlot()
# Plot(baseObj, axes = "L", ylab = "Something meaningful") 
Plot(baseObj, baseObj2, addTo = "scatterplot", axes = TRUE) 
newPoints <- rnorm(10)
newPoints2 <- newPoints * 1.2 + rnorm(10)
Plot(newPoints, newPoints2, addTo = "scatterplot", col = "red")

  # can change colour palette
  clearPlot()
  Plot(landscape) # original
  
  # can use RColorBrewer if installed
  habQualCols <- if (requireNamespace("RColorBrewer")) {
    RColorBrewer::brewer.pal(9, "Spectral")
  } else {
    colorRampPalette(c("blue", "purple"))(9)
  }
  
  mapColours <- list(
    DEM = topo.colors(50),
    forestCover = colorRampPalette(c("blue", "orange", "purple", "red"))(50),
    habitatQuality = habQualCols
  )
  setColors(landscape, n = 50) <- mapColours
  Plot(landscape, new = TRUE) # oh, how pretty!

clearPlot()
Plot(caribou, gpAxis = gpar(cex = 0.4), size = 1)

Plot(maps$DEM, gpText = gpar(cex = 0.4))

clearPlot()
Plot(maps$DEM, caribou, gpText = list(gpar(cex = 2), gpar(cex = 0.1)), new = TRUE)

# x axis gets cut off in pdf and html
clearPlot()
Plot(maps$DEM)

clearPlot()
Plot(maps$DEM, visualSqueeze = 0.6, new = TRUE)

clearPlot()
Plot(maps$DEM, legendRange = c(0, 500), new = TRUE)

Plot(maps$DEM, zoomExtent = terra::ext(c(-1, 10, -1, 20)), new = TRUE)

clearPlot()
Plot(maps$DEM)
Plot(sl, addTo = "maps$DEM", length = 0.2)

3 Modularity

One of the main purposes of the Plot function is modularity. The goal is to enable any SpaDES module to be able to add a plot to the plotting device, without being aware of what is already in the plotting device. To do this, each open plotting device has a corresponding .quickPlot object in the .quickPlotEnv, that gets created when a first Plot function is called. This object keeps the layer names, their extents, and whether they were in a RasterStack (and a few other things). So, when a new Plot is called, and new is used, then it will simply add the new layer. There may not be space on the plot device for this, in which case, everything will be replotted in a new arrangement, but taking the original R objects. This is different than the grid package engine for replotting. That engine was not designed for large numbers of plots to be added to a region; it slows down immensely as the number of plots increases.

clearPlot()
Plot(landscape)
Plot(landscape, new = list(TRUE, FALSE, TRUE), axes = "L")

clearPlot()
Plot(landscape)

# can add a new plot to the plotting window
Plot(caribou, new = FALSE, axes = FALSE)

clearPlot()
Plot(landscape)

# can add a new plot to the plotting window
Plot(caribou, new = FALSE, axes = FALSE)

clearPlot()
Plot(landscape)
landscape$forestCover[] = ((landscape$forestCover[] + 10) %% 30)
# can add a new plot to the plotting window
Plot(landscape, new = FALSE)
# note: zeros are treated as no colour by default.
# if this is not the correct behaviour, can use `zero.color`
Plot(landscape, new = FALSE, zero.color = "blue")

4 Plotting Speed

A second main purpose of the Plot function is to plot as fast as possible so that visual updates, which may be frequent, take as little time as possible. To do this, several automatic calculations are made upon a call to Plot. First, the number of plots is compared to the physical size of the device window. If the layers are RasterLayers, then they are subsampled before plotting, automatically scaled to the number of pixels that would be discernible by the human eye. If the layer is a SpatialPoints* object, then a maximum of 10,000 points will be plotted. These defaults can be adjusted by using the speedup argument. Broadly, speedup is a positive real number (default is 1). Values greater than 1 should plot faster; numbers between 0 and 1 should plot slower. There will be exceptions to this (see next), so trial and error for a particular situation is the only way to determine for certain what speedups will be gained.

system.time(Plot(landscape, caribou, maps$DEM, new = TRUE))

##    user  system elapsed 
##   0.045   0.001   0.046

system.time(Plot(landscape, caribou, maps$DEM, speedup = 10, new = TRUE))

##    user  system elapsed 
##   0.091   0.001   0.092

# can add a new plot to the plotting window

5 Overplotting: addTo

There are times when it is useful to add a plot to a different plot with a different name. In these cases, the new argument will not work. The argument addTo will allow plotting of a RasterLayer or SpatialPoints* object on top of a RasterLayer, that does not share the same name. This can be useful to see where agents are on a RasterLayer, or if there is transparency on a second RasterLayer, it could be plotted on top of a first RasterLayer.

clearPlot()
Plot(landscape)
Plot(caribou, addTo = "landscape$DEM", size = 2)

6 Clearing Plot

A device and all the meta-data associated with that device can be cleared with clearPlot function. Sometimes, the plot state can be “illegal” because it was interrupted during a previous plotting attempt. This function can be used to clear the device and all the associated metadata manually. Use “force = TRUE” to close the device, and reopen it, which is a sure-fire way to clear everything.

clearPlot()
Plot(caribou)

7 Plotting in RStudio

The built-in RStudio plot window is particularly slow. It is recommended to always create a new plotting device whenever real simulations are being done and they will be substantially faster.

Until this is fixed in RStudio, we have created a function, dev(x) which will add devices up to the number in the parenthesis, or switch to that device if it is already open. If an RStudio plot has not been called, dev() will create a new device outside RStudio. dev() can also be used with a number, as in dev(x=4) and this function will activate that device, if already exists, or open up all devices required in sequence from current device to x. Since RStudio uses 2 and 3 (for the png), a call to dev(4) will either create 3 new devices (device numbers 2, 3 and 4 because device number 1 is never used in R), or 1 new device if plotting to RStudio has already occurred.

To set the default graphics device to avoid using the RStudio one, use dev.useRSGD(FALSE), which sets the device option for your platform in order to bypass the RStudio graphics device. To reset this back to the RStudio device, use dev.useRSGD(TRUE).

8 Interacting with plots

clearPlot()
Plot(landscape)
clickValues(3) # click at three locations on the Plot device

clearPlot()
Plot(landscape)
clickExtent() # click at two locations on the Plot device

9 Replotting or redrawing the entire plot window

rePlot can be used to redraw the entire plot window, either in the same frame as before (the default) or in a new frame. If in a new frame, the new device will be created, then the entire Plot window will be redrawn in the new device. This function can also be used to redraw everything after a window gets resized manually and objects become distorted. Some plot arguments can be passed via rePlot, such as visualSqueeze, axes, but not all. Further functionality like this will likely be in a future update.K

rePlot()
rePlot(4)
rePlot(visualSqueeze = 1, axes = FALSE)
rePlot(visualSqueeze = 0.7, axes = FALSE, cols = "Reds")

10 Assigning Plot() to object

The output of a Plot call can be assigned to an object for future plotting. Simply pass the object name as if it were a Plot object. Any data changes made to the objects contained within the plot, if they are in the .GlobalEnv or a envir will be reflected in the new plot, because the data are replotted from their original source. This provides the same behaviour as a call to rePlot().

clearPlot()
plots <- Plot(landscape)

# change values
landscape$forestCover[landscape$habitatQuality > 0.9] <- 0

Plot(plots)
# same as:
rePlot()

# but can be combined with other objects
Plot(caribou, plots)
Plot(caribou, plots, sl)