Statistical Data Visualization: Raster Data

Reading, Recording, Rmarkdown

library("raster")
library("RStoolbox")
library("ggplot2")
library("dplyr")
theme_set(theme_bw())

The raster data format is used to store spatial data that lie along regular grids. The values along the grid are stored as entries in the matrix. The raster object contains metadata that associates each entry in the matrix with a geographic coordinate.
Since the data they must lie along a regular grid, rasters are most often used for continuously measured data, like elevation, temperature, population density, or landcover class.
We can create a single layer raster using the raster command. The code block below loads an elevation map measured by the space shuttle.

f <- system.file("raster/srtm.tif", package = "spDataLarge")
zion <- raster(f)

Typing the name of the object shows the metadata associated with it (but not the actual grid values). We can see that the grid has 457 rows and 465 columns. We also see its spatial extent: The minimum and maximum longitude are both close to -113 and the latitudes are between 37.1 and 37.5. A quick google map search shows that this is located in Zion national park.

zion

class      : RasterLayer 
dimensions : 457, 465, 212505  (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333  (x, y)
extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/srtm.tif 
names      : srtm 
values     : 1024, 2892  (min, max)

plot(zion)

The raster command also lets us create raster objects from scratch. For example, the code below makes a raster with increasing values in a 6 x 6 grid. Notice that we had to give a fake spatial extent.

test <- raster(
  nrows = 6, ncols = 6, res = 0.5, 
  xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
  vals = 1:36
)

plot(test)

Real-world rasters typically have more than one layer of data. For example, you might measure both elevation and slope along the same spatial grid, which would lead to a 2 layer raster. Or, for satellite images, you might measure light at multiple wavelengths (usual RGB, plus infrared or thermal for example).
Multi-layer raster data can be read in using brick. You can refer to particular layers in a multi-layer raster by using subset.

f <- system.file("raster/landsat.tif", package = "spDataLarge")
satellite <- brick(f)
satellite

class      : RasterBrick 
dimensions : 1428, 1128, 1610784, 4  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs 
source     : /Library/Frameworks/R.framework/Versions/4.0/Resources/library/spDataLarge/raster/landsat.tif 
names      : landsat.1, landsat.2, landsat.3, landsat.4 
min values :      7550,      6404,      5678,      5252 
max values :     19071,     22051,     25780,     31961

Base R’s plot function supports plotting one layer of a raster at a time. To plot more than one layer in a multichannel image (like ordinary RGB images) you can use the plotRGB function.

plotRGB(satellite)

Sometimes, it’s useful to overlay several visual marks on top of a raster image. In this case, the ggRGB function from the RStoolbox package can be used.

satellite <- projectRaster(satellite, crs = 4326)
ggRGB(satellite) +
  geom_point(
    data = data.frame(x = -113, y = 37.3), 
    aes(x = x, y = y), 
    col = "red", size = 10
  )

If we want to visualize just a single layer using ggplot2, we can use geom_raster. This assumes that the data are in “tall” format, which can be accomplished by (1) subsetting to the relevant layer and (2) converting it to a data.frame. We use xy = TRUE to remember the original locations of each pixel.

sat_df <- subset(satellite, 1) %>%
  as.data.frame(xy = TRUE)
ggplot(sat_df) +
  geom_raster(aes(x = x, y = y, fill = landsat.1)) +
  scale_fill_gradient(low = "white", high = "black") +
  coord_fixed()