I recently released an update to voluModel on CRAN. We’re now up to version 0.2.2. Some of the stuff is pretty boring but necessary (fixing broken URLs, removing a few last vestiges of the sf package, RIP). What I want to highlight is a new function that plots a vertical transect across a 3D raster map (like looking through a fish tank, as co-author Carsten Rahbek says).

Let’s get started.

Setup

For this tutorial, I will use two packages: voluModel (of course); terra for a little additional help with spatial data; and ggplot2, ggpubr, ggplotify, cowplot, and gridExtra for plotting sugar.

# August 2024 voluModel v0.2.2 mini-tutorial
# R version version 4.4.1 (2024-06-14) -- "Race for Your Life"
# Copyright (C) 2024 The R Foundation for Statistical Computing
# Platform: aarch64-apple-darwin20 (64-bit)

#Load the libraries you'll need.
library(voluModel)
library(terra)
library(ggplot2)
library(ggpubr)
library(ggplotify)
library(cowplot)
library(gridExtra)

Vertical visualization

Visualizing 3D diversity estimates (and all kinds of other 3D-structured geographic data) is a challenge. The transectPlot() was designed to help. Imagine looking at a layer cake from the top down, then slicing through it and looking at it from the side. Except in this case, the cake is a depth-structured stack of spatRaster layers.

Let’s start by loading a multi-layer raster. For simplicity, I will be simulating a datset using the smoothed dissolved oxygen dataset (see the voluModel raster processing tutorial) that comes packaged with voluModel. Note that any stacked raster can be used this way, as long as each layer is named with the depth it represents. The function was originally designed with 3D species distributions and richness estimates in mind, and may behave unpredictably when confronted with negative numbers. I welcome suggestions for how to fix this–check out the .

# Get raster
dist <- rast(system.file("extdata/oxygenSmooth.tif", 
                                 package='voluModel'))
dist <- classify(dist, cbind(-Inf, 0, 0))

# How are the layers structured?
paste0("This dataset contains ", nlyr(dist), 
       " layers, at the following depths: ", paste(names(dist), collapse = ", "))
## [1] "This dataset contains 102 layers, at the following depths: 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, 350, 375, 400, 425, 450, 475, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200, 1250, 1300, 1350, 1400, 1450, 1500, 1550, 1600, 1650, 1700, 1750, 1800, 1850, 1900, 1950, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4100, 4200, 4300, 4400, 4500, 4600, 4700, 4800, 4900, 5000, 5100, 5200, 5300, 5400, 5500"
# Show what the first 6 layers look like
plot(dist[[1:6]])

Of course, you could leave it at this–a series of horizontal maps that stack on top of each other, but this dataset has 102 layers, which would be a BIG plot. What if you are interested in the vertical structure of your data? It turns out that with voluModel::transectPlot(), this is easy. This example shows what the vertical structure of the data are along a longitudinal slice (argument sampleAxis) at -70\(^\circ\) (a.k.a. 70\(^\circ\)W; argument axisValue). Note that the plot interpolates between depth slices.

transectPlot(dist, sampleAxis = "lon", axisValue = -70)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

Longitude is plotted from South to North along the x axis, and water depth (as derived from layer names, assuming they are in meters) is plotted from deep to shallow along the y axis. The gray areas are no data values.The big gray area on the left is South America, and on the right is North America. Hispaniola is the gray stripe in the middle. By default, the function only plots to the maximum depth where there are data.

Of course, you can also slice the spatRaster along latitude (this time at 10\(^\circ\)). In this example, I have also limited the slice to between -100\(^\circ\) and -70\(^\circ\) (using the transRange argument).

transectPlot(dist, sampleAxis = "lat", axisValue = 25, transRange = c(-100,-70))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

Now it’s easier to see the depth gradient along the transect. Mexico is on the left, then Gulf of Mexico with a gradual slope up to Florida, a steep dropoff in the Straights of Florida on the Atlantic side, the Bahamas, and then the open Atlantic Ocean.

transectPlot uses viridis palettes () to color the plots. Perhaps you would prefer the mako palette, and you want to specify the scale range, so it has the same scale as another plot? Also note that transectPlot() results are a ggplot, and can be adjusted accordingly. Below, I add a title.

transectPlot(dist, sampleAxis = "lat", axisValue = 25, transRange = c(-100,-70),
             option = "mako", scaleRange = c(0,400)) +
  ggtitle("Data on a Latitudinal Transect")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

One last hint

In the above examples, it may have been challenging to visualize the transect plot in its geographic context. Here’s a little bonus code showing how to do a side-by-side plot with a map of surface conditions on the left with a red line showing the location of the transect plot shown at the right. The color scale for both plots is the same.

And now I place the plots side by side for clarity.

# Putting the plots side by side
grid.arrange(horizontalPlot,
               as_grob(plottedTransect),
               nrow = 1,heights = 5)

References

Hijmans R (2024). terra: Spatial Data Analysis. R package version 1.7-78, https://CRAN.R-project.org/package=terra.

Owens HL, Rahbek, C. (2023). voluModel: Modelling species distributions in three-dimensional space. Methods in Ecology and Evolution, 14, 841–847. https://doi.org/10.1111/2041-210X.14064.

Owens HL, Rahbek C (2024). voluModel: Modeling Species Distributions in Three Dimensions. doi:10.5281/zenodo.7372599 https://doi.org/10.5281/zenodo.7372599, R package version 0.2.2, https://CRAN.R-project.org/package=voluModel.

H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.